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Abstract 

This  thesis  compared  the  Response  Surface  Methodology  (RSM)  approach  to  the  one- 
factor-at-a-time  approach  for  calibrating  the  Bioplume  II  finite-difference  simulation  model 
of  groundwater  flow,  contaminant  transport  and  biodegradation.  The  MADE-2  data  set  of 
hydrocarbon  injection  into  pristine  groundwater  at  Columbus  Air  Force  Base,  Mississippi  was 
used  in  this  research.  Because  the  simulation  includes  both  groundwater  flow  and  contaminant 
transport,  each  calibration  included  both  phases.  The  one-factor-at-a-time  approach  reduced 
the  root-mean-squared  (RMS)  error  criterion  for  the  flow  to  0.921225  feet  in  a  total  of  36  runs 
of  Bioplume.  The  RSM  approach  reduced  the  error  criterion  to  0.918875  feet  in  a  total  of  47 
runs.  The  one-factor-at-a-time  approach  was  unable  to  reduce  the  RMS  error  criterion  for  the 
transport  calibration  below  an  initial  value  of  67.1831  parts  per  billion  (ppb)  benzene  after  21 
runs  which  spanned  the  feasible  range  of  each  of  the  parameters.  The  RSM  approach  was  able 
to  reduce  the  response  to  67.0327  ppb  after  47  runs  of  Bioplume.  The  RSM  approach  allows 
the  modeler  to  identify  parametric  regions  of  improved  response  in  a  systematic  way  that 
would  be  extremely  difficult  to  find  using  the  one-factor-at-a-time  approach.  For  this  reason 
it  may  be  very  useful  for  calibration  of  Bioplume  models  to  be  used  for  research  or  long 
term  monitoring  of  a  contaminated  site,  where  extra  prediction  accuracy  may  be  needed.  The 
major  limitations  of  this  work  were  the  use  of  inefficient  full  factorial  designs  for  the  Response 
Surface  Methodology  approach  and  the  limited  improvement  possible  on  the  response  surfaces 
possibly  due  to  the  assumption  of  homogeneous  parameter  values. 
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A  COMPARISON  OF  RESPONSE  SURFACE  METHODOLOGY  AND  A 
ONE-FACTOR- AT-A-TIME  APPROACH  AS  CALIBRATION 
TECHNIQUES  FOR  THE  BIOPLUME-II  SIMULATION  MODEL  OF 
CONTAMINANT  BIODEGRADATION 


L  Introduction 

L 1  Motivation  for  Research 

Groundwater  simulation  models  are  currently  used  to  predict  the  effectiveness  of  various 
remediation  schemes  or  the  potential  environmental  and  health  threats  posed  by  plumes  of 
groundwater  contamination.  These  models  require  parameter  estimates  and  initial  conditions 
in  the  form  of  field  data  before  they  can  predict  the  future  state  vectors  of  groundwater 
hydrology  and  contaminant  distribution.  Groundwater  models  run  with  poor  estimates  of 
parameter  values  may  result  in  inaccurate  predictions  of  aquifer  behavior,  which  may  lead  to 
counterproductive  management  decisions,  wasted  time  and  money,  and  even  environmental 
and  health  consequences. 

These  problems  might  be  avoided  if  the  models  used  were  better  calibrated.  Groundwa¬ 
ter  model  calibration  is  a  procedure  for  selecting  parameter  values  which  result  in  acceptable 
prediction  accuracy.  Professionals  might  be  more  willing  to  spend  time  on  calibration  if 
systematic  and  effective  procedures  were  outlined  for  model  calibration. 

Response  Surface  Methodology  (RSM)  has  been  proposed  as  a  calibration  technique 
for  groundwater  contamination  models  (1)  (19).  RSM  has  not  been  compared  to  any  other 
technique  to  determine  whether  it  may  represent  a  real  improvement  in  terms  of  how  many 
runs  might  be  required  for  the  method  to  converge  or  how  closely  it  may  calibrate  a  model. 
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1.2  Research  Objectives 

The  goal  of  this  research  was  to  compare  the  Response  Surface  Methodology  approach 
to  the  one-factor-at-a-time  approach  for  calibrating  the  Bioplume  II  simulation  model.  This 
research  was  accomplished  using  the  M  ADE-2  data  from  a  natural  gradient  tracer  experiment  at 
Columbus  Air  Force  Base,  Mississippi  (51).  The  measures  used  to  compare  the  two  calibration 
techniques  were  ease  of  use,  number  of  runs  required  until  improvement  in  response  slowed 
significantly,  and  the  best  response  found  using  each  method. 
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IL  Literature  Review:  The  Inverse  Problem 


In  general,  simulation  modeling  involves  the  prediction  of  outcomes  from  a  system. 
Many  models  may  be  formulated  mathematically  as: 

Z,  =  /(Zo,fi) 

where  X_^  is  the  future  state  vector  of  the  system,  is  a  vector  representing  the  current  state 
of  the  system  and  ^  is  a  vector  representing  system  parameters.  How  accurately  the  model 
predicts  future  states  of  the  system  is  a  function  of  the  parameter  set  chosen  to  run  the  model. 
In  groundwater  modeling,  as  in  many  disciplines,  direct  measurement  of  parameter  values 
may  be  difficult,  expensive,  or  even  impossible. 

The  problem  of  parameter  estimation  is  the  inverse  problem  of  mathematics,  where  the 
parameter  vector  is  solved  for  as  the  unknown.  The  inverse  problem  is  often  ill-posed  in  the 
sense  that  it  does  not  lead  to  unique  solutions  and  is  unstable  because  small  measurement 
errors  result  in  large  errors  in  parameter  estimates  (63:95).  The  nonuniqueness  of  the  solution 
suggests  that  parameter  sets  that  do  not  represent  actual  environmental  parameter  values  may 
lead  to  successful  calibration  and  even  predict  system  response. 

Anderson  and  Woessner  (2:8)  define  calibration  as  the  process  of  adjusting  parameter 
values  until  model  predictions  match  field  data.  This  process  involves  minimizing  the  differ¬ 
ence  between  the  predicted  final  state  of  the  system  and  the  actual  final  state.  This  process 
requires  the  use  of  a  historical  data  set,  i.e.  one  which  includes  initial  and  final  conditions 
for  the  site  under  consideration.  Future  states  of  the  same  system  may  be  predicted  after 
calibration  has  been  completed.  For  example  if  data  exists  for  years  one  and  two,  then  the 
model  may  be  calibrated  using  year  one  data  as  the  initial  condition  and  year  two  data  as 
the  final  condition.  Modeled  values  of  year  two  data  could  then  be  compared  to  actual  year 
two  data  and  the  difference  between  the  two  reduced  via  alteration  of  the  input  parameters. 
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Once  this  error  has  been  minimized,  the  state  of  the  system  at  year  three  may  be  predicted.  In 
groundwater  modeling,  calibration  is  also  called  history  matching. 

A  variety  of  solution  techniques  have  been  developed,  including  trial  and  error,  direct 
methods,  indirect  methods,  and  the  geostatistical  approach. 

2. 1  Trial  and  Error 

The  trial  and  error  approach  is  the  most  common  technique  in  practice  (16:199).  This 
technique  starts  with  an  approximation  of  the  parameter  values  as  either  fixed  or  spatially  vary¬ 
ing.  The  approximations  may  come  from  field  measurements,  literature  values  or  experience. 
Bair  and  others  (7:887-888)  assumed  literature  values  of  effective  porosity  and  then  used  trial 
and  error  adjustments  for  hydraulic  conductivity.  Rifai  and  others  (47:1021-1023)  performed 
calibration  on  the  Bioplume  II  groundwater  simulation  model  assuming  constant  hydraulic 
conductivity  and  recharge,  obtaining  aquifer  thickness  from  well  logs,  and  estimating  the 
reaeration  from  the  literature. 

Trial  and  error  is  probably  the  dominant  method  due  to  its  mathematical  simplicity,  but 
it  may  not  be  the  best  method  overall.  Keidser  and  Rosjberg  (32:2219)  note  that  the  skill  of 
the  modeler  plays  a  major  role  in  efficiency  and  effectiveness  of  this  method.  Carrera  and 
Neuman  (16:199)  note  the  method  is  time  consuming  and  subjective. 

2.2  Direct  Methods 

Yeh  (63:96)  defines  the  direct  method  (as  classified  by  Neuman)  or  the  equation  error 
criterion  (as  classified  by  Chavent)  to  be  the  direct  solution  of  flow  and  transport  equations  for 
parameters  using  measured  and/or  estimated  data  such  as  heads  and  concentrations.  Willis  and 
Yeh  (57:352)  state  that  this  approach  involves  an  explicit  solution  to  an  inverse  boundary  value 
problem  where  the  parameter  vector  is  treated  as  a  dependent  variable.  Wise  and  Charbeneau 
(60:429)  demonstrate  the  use  of  a  direct  semianalytical  method  to  calculate  parameter  values 
based  on  field  data.  Butcher  and  Gauthier  (14:73)  applied  a  direct  analytical  solution  to  the 
inverse  boundary  value  problem  of  determining  the  source  geometry  of  a  non-aqueous  phase 
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layer.  Statistical  methods  and  iterative  approaches  may  be  coupled  with  direct  methods,  but 
in  general  the  direct  methods  are  characterized  by  their  explicit  solutions  and  lack  of  iterative 
technique.  Their  disadvantage  is  that  the  cumbersome  mathematics  behind  the  direct  methods 
makes  them  awkward  for  solving  complex  problems. 


2.3  Indirect  Methods 


Indirect  methods,  also  called  Output  Error  Criterion  methods  by  Chavent  refer  to  an 
iterative  minimization  of  a  norm  of  observed  versus  calculated  data  in  order  to  approximate 
model  parameters  (63:96).  Trial  and  error  may  be  classified  as  an  indirect  method  because 
it  involves  this  successive  approximation  approach.  Most  indirect  methods  have  an  error 
criterion  by  which  the  goodness-of-fit  of  the  parameter  estimation  is  judged.  Anderson  and 
Woessner  (2:238-241)  describe  three  calibration  criteria  for  matching  simulated  to  observed 
groundwater  head  values:  Mean  Error  (ME),  Mean  Absolute  Error  (MAE),  and  Root  Mean 
Squared  Error  (RMS).  The  Mean  Error  is  the  average  difference  between  simulated  and 
observed  values  at  data  points  (well  nodes).  Mean  Absolute  Error  is  average  absolute  value 
of  the  difference  between  simulated  and  observed  values.  Root  Mean  Squared  Error  is  the 
square  root  of  the  average  squared  difference  between  simulated  and  observed  values  of  head. 
The  criteria  are  as  shown  in  equations  1,  2,  and  3  where  n  is  the  number  of  observation  or 
well  nodes,  and  the  subscripts  o  and  s  represent  observed  and  simulated  values,  respectively. 


ME 

MAE 

RMS 


^  f^iho  -  K)i 


”7=1 


-f2\iho-hs)i 


n 


i=l 


ztiK-Kn 


”7=1 


(1) 

(2) 

(3) 


A  response  surface  can  be  formed  by  the  functional  relationship  between  the  parameter 
values  and  the  error  criterion  selected.  A  relationship  of  the  form  Y  =  f{6)  would  then 
be  developed,  where  0  is  a  vector  of  parameter  values  and  Y  is  the  response  function. 
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Identification  of  an  acceptable  parameter  vector  is  performed  by  moving  along  the  response 
surface  looking  for  an  optimal  response  value. 

A  variety  of  methods  have  been  developed  to  identify  local  minima.  Devlin  (25:327) 
demonstrated  that  Simplex  optimization  may  be  used  to  identify  local  mimina  of  a  response 
surface  in  order  to  estimate  parameter  values.  Yeh  and  others  (62:35)  developed  an  iterative 
(indirect)  parameter  estimation  routine  involving  two  steps:  ( 1 )  cokriging  to  estimate  unknown 
parameter  values  and  (2)  steady-state  simulation  modeling  with  observed  heads  as  constant, 
boundary  values.  Doughty  and  others  (26:1741)  developed  an  indirect  technique  to  solve 
the  inverse  problem  based  on  the  use  of  fractals.  Carrera  and  Neuman  (17:218)  describe  a 
method  to  calculate  optimal  gradients  to  employ  in  the  iterative  method  using  a  finite  element 
approximation  of  the  gradient  of  a  least  squares  criterion  including  a  head  residual  and  a 
penalty  function.  Olsthoom  (42:44)  proposed  an  indirect  method  to  estimate  parameter  values 
by  generating  log  parameter  vectors,  running  the  model,  computing  the  error  criterion,  a 
Jacobian  matrix,  and  a  gradient  vector  of  the  error  criterion.  The  algorithm  iteratively  updates 
the  parameter  vector  after  each  step.  Freeze  (27:751)  noted  that  several  statistical  methods 
may  be  used  to  support  indirect  solutions  to  the  inverse  problem  including  weighted  least 
squares  estimation,  bayesian  estimation,  and  maximum  likelihood  estimation.  Adams  (1)  and 
Cottman  (19)  demonstrated  that  Response  Surface  Methodology  may  be  used  to  solve  the 
inverse  problem  of  flow  model  calibration  by  locating  minima  on  a  response  surface  using  the 
sum  of  squared  error  criterion. 

2.4  Geostatistical  Approach 

The  geostatistical  approach  involves  kriging  and  cokriging  to  estimate  parameter  values. 
Kriging  is  a  statistical  technique  for  optimal  estimation  of  a  spatially  distributed  parameter  for 
which  some  values  have  been  measured  at  a  discrete  points  (observation  wells).  Cokriging 
is  a  related  technique  by  which  an  unmeasured  parameter’s  values  are  estimated  from  the 
values  of  a  measured  parameter  which  is  correlated  to  the  parameter  to  be  estimated  (24:323). 
Either  of  these  techniques  may  be  used  to  estimate  parameter  values  from  limited  field  data. 
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These  techniques  may  also  be  used  to  fill  data  gaps  where  needed  in  connection  with  another 
calibration  technique.  Sun  and  others  (52:90)  used  kriging  to  estimate  parameters  as  contin¬ 
uously  varying  via  measured  point  values  and  the  use  of  a  covariance  matrix.  Keidser  and 
Rosjberg  (32:2230-2231)  compared  four  approaches  to  solving  the  inverse  problem:  kriging 
with  prior  data  on  response,  kriging  with  zonation  (establishing  areas  of  constant  parameter 
value),  kriging  alone,  and  zonation  alone.  They  found  that  kriging  with  zonation  resulted 
in  the  best  model  fit,  except  when  data  was  scarce  or  there  was  large  measurement  error,  in 
which  case  zonation  alone  was  best.  Freeze  and  others  (27:751)  explained  that  geostatistical 
approaches  differ  from  other  statistical  methods  to  solve  the  inverse  problem  because  they 
view  parameters  as  spatially  stochastic.  Carrera  and  Glosorio  (15:281)  concluded  that  indirect 
statistical  methods  are  superior  to  geostatistical  methods. 

2.5  Sensitivity  Analysis 

As  stated  earlier,  the  inverse  problem  is  often  unstable  because  small  measurement  errors 
generally  result  in  large  errors  in  parameter  estimates  (63:95).  Errors  in  parameter  estimates 
generally  produce  errors  in  predicted  system  response.  Sensitivity  analysis  of  parameter 
estimates  is  essential  to  determine  the  effect  of  parameter  uncertainty  on  model  prediction. 

Sensitivity  coefficients  are  the  partial  derivatives  of  response  (head  error  or  concentra¬ 
tion  error)  with  respect  to  each  of  the  model  parameters.  Willis  and  Yeh  (57:376)  describe 
three  methods  of  computing  sensitivity  coefficients:  the  influence  coefficient  method,  the  sen¬ 
sitivity  equation  method,  and  the  variational  method.  The  influence  coefficient  method  was 
developed  by  Becker  and  Yeh  (57:376)  and  involves  altering  parameter  values  and  succes¬ 
sively  resolving  the  simulation  model  to  observe  the  effect  on  model  response.  The  sensitivity 
equation  method  involves  taking  the  partial  derivatives  of  the  flow  and/or  transport  equations 
with  respect  to  each  of  the  parameters  and  numerically  solving  for  the  sensitivity  coefficients 
(57:377).  The  variational  method,  developed  by  Jacquard  and  Jain,  Carter  and  others,  and 
Sun  and  Yeh  involves  numerically  estimating  the  coefficients  using  a  triple  integral  over  the 
subdomain  of  each  node  (57:378). 
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Other  formal  approaches  have  been  developed  to  analyze  sensitivity  without  the  use 
of  sensitivity  coefficients.  Keidser  and  Rosjberg  (32:2219)  suggest  that  if  model  parameters 
are  random  variables  then  statistical  estimation  and  inference  may  be  performed  to  assess 
reliability.  Brooks  and  others  (13:2996-2997)  proposed  a  sensitivity  analysis  process  for 
groundwater  modeling  based  on  using  a  simplex  approach  to  identify  extreme  cases  of  model 
response  based  on  variations  in  parameter  values  between  feasible  limits.  Spear  and  others 
(49:3161)  developed  a  method  of  parametric  sensitivity  analysis  involving  monte  carlo  simu¬ 
lation  that  allows  variation  of  more  than  one  parameter  at  a  time.  Yeh  (63:103)  cited  Yeh  and 
Yoon,  and  Shah  and  others  who  noted  parameter  uncertainty  may  be  measured  as  the  norm  of 
the  covariance  matrix  of  estimated  parameters. 

2.6  Parameterization 

Parameterization  refers  to  the  discretization  of  actual  environmental  variables  into  a 
given  number  of  model  parameters.  The  parallel  concept  of  zonation  is  the  assignment 
of  zones  of  constant  parameter  value  to  model  parameters.  Carrera  and  Neuman  (16:206) 
note  that  too  many  zones  may  increase  error  in  parameter  estimates.  As  parameterization 
increases  system  modeling  error  decreases,  but  parameter  uncertainty  increases  (63:95).  Sun 
and  others  (52:89)  noted  that  because  overparameterization  may  increase  the  variance  of 
identified  parameters,  optimal  parameterization  is  a  balance  between  the  error  criterion  (such 
as  the  residual  of  least  squares  error)  and  parameter  uncertainty  (measured  for  example  by  the 
covariance  matrix  of  the  estimated  parameters). 

Several  approaches  to  parameterization  have  been  developed.  The  simplest  technique 
is  to  establish  lumped  parameter  values  which  ignore  spatial  variations  and  assume  the  entire 
aquifer  to  be  homogeneous  (44:2).  Yeh  (63:98)  explained  that  parameterization  may  be  by 
zonation,  where  parameter  values  are  constant  over  each  zone,  or  by  interpolation,  where 
parameter  values  vary  between  nodes.  In  the  finite  element  interpolation  method  (52:90) 
the  parameter  varies  continuously  via  the  basis  functions  which  apply  between  sets  of  nodal 
values.  In  the  stochastic  inverse  method  (52:90)  the  parameter  varies  over  a  field  via  statistical 
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parameters  and  is  estimated  with  the  maximum  likelihood  approach  combined  with  cokriging. 
Sun  and  others  (52:101)  propose  parameterizing  the  model  in  connection  with  the  geological 
characterization  of  the  subsurface. 

2.7  Summary 

Extensive  research  has  been  completed  on  groundwater  simulation  model  calibration. 
The  problem  is  that  many  of  these  techniques  have  not  been  evaluated  for  efficiency.  If  a 
method  is  no  more  efficient  than  the  trial  and  error  approach,  may  not  find  support  among 
modelers.  Response  Surface  Methodology  provides  an  example.  This  technique  has  been 
shown  to  work  for  calibration  of  groundwater  simulation  models  (1)  (19),  but  has  not  been 
demonstrated  for  contamination  models  or  compared  to  other  techniques.  This  thesis  will 
evaluate  the  Response  Surface  Methodology  to  examine  whether  it  may  be  more  efficient  or 
more  effective  than  the  baseline;  one-factor-at-a-time  or  trial  and  error  approach. 
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III.  Background 


3.1  Bioplume  II 

Bioplume  II  is  a  FORTRAN  simulation  model  of  the  biodegradation  of  a  contaminant  in 
groundwater.  It  was  written  by  Dr.  Hanadi  Rifai  (46)  (47)  of  Rice  University  as  an  alteration 
of  the  Method  of  Characteristics  (MOC)  finite  difference  model  of  fate  and  transport  written  by 
Konikow  and  Bredehoeft  (34)  at  the  USGS.  The  model  simulates  steady  state  or  transient  flow 
and  both  advective  and  dispersive  transport,  as  well  as  both  dissolved  oxygen  and  hydrocarbon 
plumes.  At  the  end  of  each  time  step,  Bioplume  II  simulates  the  reaction  between  the  two 
plumes  using  a  mineralization  ratio  defined  by  the  user.  This  approach  differs  from  one  where 
biodegradation  is  modeled  using  a  simple  decay  term.  The  model  is  a  standard  in  the  field 
of  environmental  restoration  and  is  commonly  used  to  design  remediation  systems  that  utilize 
natural  attenuation  or  injection  of  oxygenated  water  (48:55). 

Bioplume  II  requires  an  input  file  which  includes  all  the  data  that  the  model  needs  to 
simulate  flow  and  transport.  This  input  file  contains  values  of  constant  and  spatially  varying 
parameters,  such  as  porosity  and  hydraulic  conductivity,  as  well  as  initial  concentrations  of 
oxygen  and  hydrocarbon  contaminant.  When  Bioplume  II  is  run,  the  system  prints  changing 
values  of  plume  concentrations  and  other  data  related  to  the  modeled  groundwater  system. 

Bioplume  II  requires  two  layers  of  calibration  because  there  are  two  layers  of  simu¬ 
lation.  The  hydraulic  flow  and  the  transport  equations  must  both  be  considered.  First  a  set 
of  flow  parameters  are  adjusted  until  the  predicted  hydraulic  heads  match  the  expected  heads 
sufficiently.  Second,  a  set  of  transport  parameters  are  adjusted  until  the  predicted  concentra¬ 
tions  are  close  to  the  expected  final  concentrations.  After  both  steps  are  accomplished  the 
calibration  is  complete. 

3.2  Calibration  Techniques 

Two  calibration  techniques  were  compared  in  this  thesis:  a  one-factor-at-a-time  ap¬ 
proach  and  the  Response  Surface  Methodology  (RSM)  technique.  The  RSM  technique  in- 
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volves  several  procedures  which  can  be  used  to  locate  the  optimal  response  on  an  empirical 
response  surface.  RSM  will  be  evaluated  by  comparing  the  number  of  runs  required  for 
calibration  and  the  final  degree  of  calibration  via  the  Root  mean  Squared  error  criterion. 
The  one-factor-at-a-time  approach  involves  improving  on  the  response,  or  fit,  by  iteratively 
changing  one  parameter  at  a  time. 

3.2.1  One  Factor  at  a  Time  Calibration.  The  trial  and  error  approach  starts  with 
an  educated  guess  at  each  of  the  parameter  values  and  then  the  modeler  manually  alters  the 
parameters,  one  at  a  time,  until  an  acceptable  match,  or  fit,  between  the  predicted  and  actual 
final  conditions  is  achieved.  In  practice,  trial  and  error  is  not  any  one  standard  approach,  but 
the  set  of  all  manual  iterative  guessing  calibration  techniques.  For  the  purposes  of  this  work 
the  trial  and  error  technique  means  the  one-factor-at-a-time  approach. 

To  implement  this  approach,  the  parameters  are  listed  in  the  order  in  which  they  are  to 
be  improved  upon.  Next  the  Bioplume  II  model  is  run  at  an  initial  parameter  setting.  The 
response  for  this  parameter  set  is  the  error  criterion  for  that  point.  The  first  parameter  is  then 
increased  by  a  set  amount,  equal  to  the  first  significant  digit  (e.g.  32.7  would  be  increased  to 
42.7  or  decreased  to  22.7)  and  the  model  run  again.  If  the  initial  parameter  value  is  zero  any 
reasonable  step  size  may  be  chosen. 

If  the  response  does  not  improve,  the  parameter  is  decreased  by  the  set  amount  to  test 
in  the  other  direction.  If  no  improvement  in  response  is  found  in  either  direction,  then  the 
parameter  should  be  multiplied  and  then  divided  by  increasing  powers  of  ten  until  either  a 
change  in  response  is  seen  or  the  limit  of  a  reasonable  range  for  that  parameter  is  reached. 
When  the  parameter  improves,  it  is  altered  again  in  the  direction  of  improving  response. 

This  process  is  continued  until  response  stops  improving.  The  last  three  settings  make 
two  ranges  which  are  then  split  in  half.  For  example  if  the  model  was  run  for  parameter  x  at 
a;  =  1,  a;  =  2  and  x  =  3  and  improved  between  one  and  two,  but  error  increased  for  three, 
then  the  next  step  would  be  to  split  the  two  ranges  and  run  the  model  at  x  =  1.5  and  x  =  2.5. 
The  better  of  these  two  responses  is  chosen  for  further  exploration  with  the  ranges  being  split 
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in  half  at  each  step.  Whenever  two  responses  are  equal,  the  choice  of  which  one  to  further 
explore  should  be  made  randomly. 

Once  improvement  in  response  ceases,  the  parameter  is  set  at  the  last  value  for  which 
improvement  had  continued  and  the  procedure  is  repeated  for  the  next  parameter.  This 
procedure  is  continued  until  the  RMS  response  goes  below  some  established  cut  off  value.  If 
all  parameters  are  improved  without  meeting  a  termination  criterion,  the  modeler  goes  through 
the  list  again  beginning  with  the  first  parameter  improved  upon.  If  no  improvement  is  seen 
after  going  through  the  whole  list  of  parameters  then  the  method  may  be  seen  as  stalled. 

3.2.2  Response  Surface  Methodology.  Response  Surface  Methodology  (RSM)  is  a 
technique  for  finding  local  optimal  response  values  from  an  empirical  model.  The  response  in 
this  case  is  a  function  of  the  difference  between  the  predicted  and  actual  final  condition.  The 
actual  surface  that  RSM  explores  is  the  n-dimensional  surface  formed  by  assigning  a  response 
value  to  each  coordinate  point  (^i  ,^2,  •  •  •  ,  6n),  where  Oi  is  the  value  of  the  ith  factor.  Note 
that  the  value  of  this  response  surface  could  be  determined  at  any  point  by  running  the  model 
using  the  parameter  values  corresponding  to  the  coordinate  (^i,  ^2?  •  •  •  >  ^n)  and  calculating 
the  error  criterion  for  that  run.  This  process  allows  the  modeler  to  explore  the  surface  at  any 
point  without  knowledge  of  the  nature  of  the  response  function  Y  —  f{0)  itself.  The  RSM 
technique  allows  the  modeler  to  examine  a  series  of  local  approximations  of  the  surface  and 
move  towards  parameter  regions  of  improved  response.  The  method  is  based  on  the  statistical 
fields  of  regression,  analysis  of  variance,  and  experimental  design.  The  actual  process  of  RSM 
starts  with  a  screening  phase. 

3.2.2. 1  Screening  Phase.  The  purpose  of  this  phase  is  to  determine  which 
parameters  significantly  influence  the  response  of  the  system.  This  process  may  be  especially 
helpful  when  there  are  a  large  number  of  potential  factors.  Screening  is  generally  accomplished 
by  conducting  experiments  according  to  a  Plackett-Burman  or  2*7/  fractional  factorial  design, 
but  may  be  performed  less  efficiently  with  a  two-level  full-factorial,  or  2^,  design.  A  two- 
level  full  factorial  design,  for  k  factors  contains  every  possible  combination  of  the  high  and 
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low  settings  for  the  k  factors  (0i,  02,  •  •  • ,  ^fc)-  For  example  if  two  parameters  were  under 
consideration,  the  model  would  be  run  at  four  different  parameter  combinations  or  coordinate 
pairs  (01, 02)  because  for  /;  =  2,  2*  =  4. 

For  convenience  and  ease  of  analysis  the  parameter  values  are  linearly  transformed  into 
coded  variables  Xk'. 


Xk 


0fc  —  0fcO 

Sk 


(4) 


where  6ko  is  the  value  of  0  at  the  center  of  the  design  region  and  Sk  is  half  the  difference 
between  the  high  and  low  levels  of  9k  (the  half-range  of  the  region).  The  values  of  the  coded 
variables  are  +1  for  the  high  level  and  -1  for  the  low  level. 

The  response  values  (error  values  from  the  experimental  runs)  are  fit  to  a  regression 
equation  in  the  fc-coded  variables  with  or  without  interaction  terms.  For  a  three  factor  case 
the  regression  equation  would  be: 


Y  =  /3o  +  ^iXi  -I-  ^2X2  +  ^zXz  ^12X1X2  +  ^IzXlXs  +  P2ZX2XZ  +  ^\2ZX\X2Xz  +  ti  (5) 

Note  that  the  error  term  in  equation  5  is  not  random  as  is  often  the  case  with  regression 
equations.  In  this  work  the  error  term  is  constant  because  the  model  is  deterministic.  If  the 
same  parameter  set  is  used  as  input  for  Bioplume  II  in  two  different  runs,  the  output  will  be 
the  same.  That  is  not  to  say  there  is  no  error:  the  response  surface  is  an  error  surface,  but  the 
error  is  completely  due  to  lack-of-fit.  Therefore  there  is  no  pure  error  and  the  error  surface  is 
not  stochastic. 

The  regression  equation  of  the  error  surface  indicates  which  factors  have  the  most 
influence  on  the  quality  of  the  model  calibration.  Identification  of  these  factors  can  be 
accomplished  using  a  stepwise  regression  procedure  or  a  normal  probability  plot.  Since  the 
deterministic  nature  of  the  data  undermines  the  assumptions  required  by  the  standard  t  and  F 
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test  statistics,  care  must  be  taken  when  determining  factor  importance.  Those  factors  accepted 
as  influential  are  then  considered  in  the  first  order  design  phase. 

3.2.2.2  First-Order  Design  Phase.  In  the  first  order  design  phase,  the 
response  surface  is  assumed  to  be  locally  linear.  RSM  uses  two  level  factorial  designs  to 
approximate  the  locally  linear  nature  of  the  surface  and  a  gradient  search  to  locate  a  region 
of  improved  response.  The  process  of  generating  an  experimental  design,  fitting  a  first-order 
linear  model,  and  conducting  a  gradient  search  continues  until  either  the  response  is  sufficiently 
low  or  significant  improvement  in  response  is  not  seen  with  further  iterations. 

If  the  first-order  design  does  not  converge  on  an  acceptably  low  error  value,  then  a 
second-order  design  may  be  needed.  Refer  to  Box  and  Draper  (12),  for  example,  for  more 
information  on  the  second-order  phase. 
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IV.  Experimentation 


4.1  Description 

The  evaluation  of  RSM  as  a  calibration  technique  involved  repeated  calibration  of  the 
Bioplume  II  model  using  both  RSM  and  the  one-factor-at-a-time  baseline  approach.  The 
calibration  exercises  involved: 

1.  Developing  the  Bioplume  II  input  file  from  the  data. 

2.  Developing  Final  Condition  Grid  from  the  data. 

3.  Performing  a  one-factor-at-a-time  calibration  on  the  simulation,  documenting  the  num¬ 
ber  of  runs  required  to  complete  the  calibration,  and  recording  the  best  response  identi¬ 
fied. 

4.  Performing  an  RSM  calibration,  documenting  the  number  of  runs,  and  recording  the 
best  response  identified. 

5.  Comparing  the  one-factor-at-a-time  approach  to  RSM. 

Cottman  (19)  and  Adams  (1)  used  the  RSM  approach  to  calibrate  a  groundwater  flow 
model.  This  work  differs  from  that  of  Cottman  and  Adams  in  that  the  RSM  calibration  was 
conducted  in  two  stages  and  compared  to  a  baseline.  The  two  stage  calibration  refers  to  the 
fact  that  two  calibration  surfaces  exist:  one  for  the  flow  calibration  and  the  other  for  the 
transport  ealibration. 

The  initial  parameter  values  chosen  for  each  phase  of  the  calibration  (i.e.  flow  or 
transport)  were  selected  from  the  feasible  range  of  each  of  the  parameters  under  consideration. 
.The  parameters  under  consideration  for  the  flow  calibration  were  the  hydraulic  conductivity, 
the  porosity,  and  the  diffuse  recharge  rate.  The  parameters  adjusted  for  transport  calibration 
were  longitudinal  dispersivity,  ratio  of  longitudinal  to  transverse  dispersivity,  stoichiometric 
ratio  (for  biotic  reaction  between  oxygen  and  benzene),  retardation  factor  and  reaeration 
coefficient. 
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The  actual  calibration  using  either  technique  required  the  use  of  an  error  criterion.  In 
this  work,  the  root  mean  square  (RMS)  criterion  was  selected  because  it  was  hoped  that  it 
would  yield  a  more  linear  response  surface  than  that  posed  by  the  SSE  criterion  used  by 
Cottman  (19). 

4.1.1  Mathematical  Setting  of  the  Problem.  The  mathematical  problem  solved  in 
this  research  is  the  inverse  problem  of  mathematics,  discussed  in  chapter  2.  Specifically,  the 
parameters  for  the  Bioplume-II  model  are  estimated  using  an  empirical  response  function.  The 
error  generated  by  estimating  model  parameters  was  minimized  using  both  the  one-factor-at- 
a-time  and  RSM  approaches. 

The  Root  Mean  Squared  Error  is  a  function  of  the  difference  between  the  simulated  final 
condition  and  the  actual  final  condition  (kriged  from  real  data  in  this  case).  The  simulated  final 
condition  may  be  found  for  any  combination  of  parameter  values,  0  by  running  the  simulation 
model,  Bioplume-II.  The  domain  of  the  problem  is  the  n  1  dimensional  space  formed  by  the 
n  factors  being  calibrated  and  the  error  response.  The  boundaries  of  the  domain  are  the  range 
of  feasible  values  of  each  paramater  found  in  the  literature.  The  initial  condition  is  provided 
by  field  data  from  the  MADE-2  site.  The  boundary  conditions  of  the  simulation  model  are 
constant  hydraulic  heads  based  on  field  data:  a  first-order  boundary  condition. 

To  simplify  the  solution  homogeneity  was  assumed  for  all  factors  that  were  calibrated. 
Isotropy  was  assumed  for  all  factors  except  dispersivity,  which  was  dealt  with  in  the  calibration. 
To  the  degree  that  these  were  poor  assumptions,  additional  error  was  introduced. 

4.2  Conduct 

The  experiment  required  construction  of  an  input  file,  determination  of  a  final  condition 
and  both  flow  and  transport  calibration  using  the  one-factor-at-a-time  approach  and  the  RSM 
approach. 
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4.2.1  Input  File  Preparation.  The  preparation  of  the  input  file  involved  obtaining 
the  MADE-2  data  set,  developing  the  finite  difference  grid,  determining  an  acceptable  time 
step,  and  constructing  the  input  file. 

4.2. 1.1  Data.  The  first  step  was  to  develop  a  Bioplume  II  input  file  from 
the  data  provided  by  Dr.  Stauffer  on  the  MADE-2  site  (51)  (50).  The  data  from  the  site 
was  measured  in  three  dimensions.  Bioplume-II  simulates  only  two  dimensions,  however, 
so  some  data  had  to  be  transformed  via  vertical  averaging.  The  data  set  included  hydraulic 
conductivity,  dissolved  oxygen,  contaminant  concentrations,  and  water  surface  elevations. 
Readings  were  available  from  June  1990  through  September  1991.  The  input  file  represents 
the  condition  at  the  site  in  June  1990.  Bioplume  was  used  to  "predict"  conditions  in  September 
1991.  These  predicted  values  were  then  compared  to  the  measured  data  from  September  1991. 
The  precision  of  the  measurements  and  the  uncertainty  in  the  data  from  the  MADE-2  site  are 
not  important  in  this  research  because  the  goal  is  to  compare  the  two  methods  of  calibrating 
the  model.  If  any  conclusions  were  to  be  drawn  from  the  field  data  itself,  then  the  quality  of 
the  data  would  be  more  important. 

Only  parts  of  the  data  were  needed.  The  hydraulic  conductivity  was  altered  as  one  of 
the  calibration  parameters,  and  was  represented  as  a  single  homogeneous  value.  Because  of 
this  simplification,  a  rough  average  value  was  all  that  was  needed  to  represent  that  parameter. 
This  simplification  would  also  lead  to  an  increase  in  the  minimum  value  of  the  error  response. 
This  had  profound  implications  for  the  comparison  between  the  two  methods  because  the 
techniques  were  limited  in  how  well  they  could  calibrate  the  model.  Because  Bioplume  II 
has  the  capability  to  model  injection  of  contaminants,  initial  plume  data  was  not  needed.  The 
initial  water  surface  profile  and  dissolved  oxygen  values  (50:89-96)  were  important,  but  before 
they  could  be  added  to  the  input  file,  a  finite  difference  grid  had  to  be  developed  onto  which 
values  could  be  placed. 

4.2. 1.2  Grid  Design.  The  maximum  grid  size  that  Bioplume-II  allows  is  20 
X  30  nodes.  This  size  was  chosen  to  increase  the  accuracy  of  the  modeling;  in  general  denser 
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grids  exhibit  reduced  error.  Due  to  the  size  of  the  site,  a  10  x  10  meter  cell  seemed  appropriate, 
but  two  factors  complicated  the  grid  spacing. 

First,  the  extent  of  plume  spreading  was  not  as  large  as  the  site.  After  benzene  was 
selected  as  the  contaminant  to  be  modeled,  the  grid  size  was  reduced  in  order  to  minimize  the 
size  of  the  modeled  region,  while  still  capturing  the  entire  benzene  plume.  Those  water  table 
and  dissolved  oxygen  data  points  outside  the  region  of  the  benzene  plume’s  maximum  extent 
were  deleted  and  the  remaining  points  were  used  to  construct  the  data  set  with  the  Surfer  (31) 
software  package. 

The  second  factor  was  that  the  Surfer  software  used  to  approximate  the  values  of 
spatially  varying  groundwater  variables  (31)  can  only  view  the  nodes  as  intersections  on  a 
grid,  not  as  cells  of  the  type  used  by  a  finite  difference  model.  To  use  Surfer  the  intersections 
must  be  viewed  as  the  centers  of  the  finite  difference  nodes.  Surfer  then  calculates  the  length 
and  width  of  the  cells  by  dividing  the  total  dimension  by  the  number  of  internodal  spaces.  For 
twenty  nodes.  Surfer  uses  twenty  grid  intersections  and  has  nineteen  internodal  spaces.  For 
thirty  nodes  Surfer  has  twenty-nine  internodal  spaces. 

The  final  length  and  width  of  the  modeled  region  was  set  automatically  by  Surfer  to  the 
range  of  well  coordinates  that  were  included  after  the  data  point  reduction  described  above. 
The  total  length  was  391 .1724  feet  and  the  total  width  was  367.357  feet,  with  the  injection  well 
at  196.144  feet  by  81.1 144  feet  from  the  comer  of  the  modeled  region.  The  placement  of  the 
well  was  based  on  the  actual  location  of  injection  wells  at  the  MADE-2  site  (50).  Therefore, 
as  shown  below,  the  longitudinal  length  of  the  nodes  was  set  to  13.49  feet  and  the  transverse 
dimension  of  the  cells  to  19.33  feet. 


(367.357  meters) (3.28  ft/m) 
19  spaces 


13.49  feet 


(391.1724  meters)(3.28  ft/m) 
29  spaces 


19.33  feet 
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The  injection  well  was  then  placed  at  1 1  nodes  in  the  x-direction  !=a  11)  and  23 

nodes  in  the  y-direction  ~  7)  from  the  comer  of  the  modeled  region. 

In  order  to  ensure  stability,  the  nodal  discretization  needed  to  be  acceptably  fine.  A 
maximum  nodal  length  in  the  flow  direction  could  be  calculated  via  the  cell  Peclet  Number 
(Pec)  as: 

p  advection  vAy  vAy  Ay 

^  dispersion  Dl  oly, 

where 

Pcc  =  Cell  Peclet  Number 

V  =  Groundwater  Velocity 

Ay  =  Longitudinal  Length  of  Cell 

ai,  =  Longitudinal  Dispersivity 

Dl  =  Longitudinal  Dispersion  Coefficient 

The  cell  Peclet  number  should  be  no  more  than  two  to  ensure  numerical  stability  (61 :65). 
As  seen  in  equation  8  the  node  length  in  the  direction  of  flow  should  be  no  more  than  twice 
the  longitudinal  dispersivity. 


Pec 

< 

2 

(6) 

Ay 

— 

< 

2 

(7) 

OtL 

Ay 

< 

2q:l 

(8) 
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Young  and  Boggs  (64;  11-15)  estimated  longitudinal  dispersivity  at  the  site  to  be  between 
10  and  42  meters.  Mercado  (39)  suggested  estimating  longitudinal  dispersivity,  a^,  as: 

ocL  =  (9) 

"  avg 

where 

Pavg  =  Average  Value  of  Permeability 
dd  =  Standard  Deviation  of  Permeability 
Ld  =  Mean  travel  Distance 

Boggs  and  others  (10:3282)  present  permeability  distributions  for  four  core  samples 
from  the  MADE-2  site.  They  also  note  the  tracer  plume  traveled  about  280  meters  in  594 
days  (10:3288).  A  value  of  36.3  feet,  calculated  via  equation  9,  is  within  the  range  prescribed 
by  Young  and  Boggs  (64:1 1-15)  though  at  the  low  end  of  that  range  (10  meters  =  32.8  feet). 
That  value  is  acceptable  however  because  a  smaller  value  of  ai,  is  more  conservative  in  terms 
of  nodal  discretization  because  it  should  lead  to  a  finer  grid.  Keeping  Ay  below  twice  the 
longitudinal  dispersivity,  the  maximum  longitudinal  length  of  a  node  is  72.6  feet.  This  value 
is  safely  above  the  value  of  13.49  to  be  used  in  this  work. 

4.2. 1.3  Time  Step.  The  length  of  the  time  step  also  impacts  numerical 
stability.  A  shorter  time  step  increases  accuracy,  but  can  increase  computer  time.  A  simple 
equation  can  be  used  to  calculate  a  maximum  time  step.  For  stability  of  the  algorithm  the  value 
of  V,  the  Courant  number,  defined  hy  v  =  where  v  equals  the  velocity  of  groundwater, 
should  be  no  more  than  one  (61:65).  Solving  for  At  in  this  equation  yields: 


V 


Given  that  v  <\\ 
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when  Ay  =  10  meters; 


.  10  meters 

At  < - 

V 

Velocity  can  be  estimated  from  measured  tracer  plume  movements  of  approximately 
160  meters  in  503  days  and  280  meters  in  594  days  for  rates  of  0.318  and  0.471  meters  per 
day  (10:3287-3288). 

Therefore: 


At  < 


10  meters 
0.471  meters/day 


(10) 


At  <  21.2  days 


(11) 


The  maximum  time  step  should  therefore  be  21.2  days  in  order  to  maintain  numerical 
stability.  Bioplume  allows  the  modeler  to  set  up  multiple  pumping  periods  which  can  each 
have  different  length  time  steps.  In  order  to  model  the  injection,  an  initial  pumping  period 
was  needed  which  coincided  with  the  injection  time.  Once  the  injection  time  was  over,  the 
remaining  time  could  be  divided  evenly  into  any  convenient  units.  Since  the  injection  took 
two  days,  two  one-day  time  steps  were  included  in  the  first  pumping  period.  The  remaining 
time  was  divided  into  five  pumping  periods  of  92  one-day  time  steps  each.  These  time  steps 
directed  the  computer  to  model  the  fifteen  month  period  with  time  steps  much  smaller  than 
the  21.2  day  maximum  required  for  numerical  stability. 


4.2. 1.4  File  Assembly.  Surfer  was  used  to  set  up  the  Bioplume-II  input 
file  by  estimating  initial  values  of  dissolved  oxygen  and  water  table  elevation  for  each  finite 
difference  node.  Other  characteristics  of  the  initial  site  condition  were  set  to  specific  values. 
The  transmissivity  was  set  to  0.00226042  square  feet  per  second  based  on  an  approximation 
of  the  hydraulic  conductivity  from  the  data  divided  by  1.96  feet  of  depth.  The  1.96  feet  was 
both  the  aquifer  screening  length  and  the  modeled  depth  of  the  aquifer.  Injection  wells  were 
set  to  values  given  in  Boggs  (9:2-1).  In  most  cases,  the  value  for  leakance,  recharge  of  water 
to  the  aquifer,  was  set  to  zero.  At  the  upstream  and  downstream  ends,  however,  the  leakance 
parameter  was  set  to  keep  the  heads  constant  for  a  first  order  boundary  condition. 

The  hydraulic  flow  modeling  was  assumed  to  be  steady  state.  To  test  this  assumption. 
Surfer  was  used  to  draw  contour  maps  of  the  water  table  surface  over  the  time  period  to  be 
modeled.  The  contour  maps  confirmed  the  acceptability  of  the  steady  state  assumption.  They 
indicated  that  for  the  time  period  under  consideration  at  the  site,  the  hydraulic  gradients  did 
not  change  significantly.  For  more  information  on  the  structure  of  the  input  file,  refer  to  Rifai 
(46).  The  initial  input  file  is  included  as  Appendix  B. 

4.2.2  Final  Condition.  To  determine  if  the  model  was  calibrated  after  a  given  run, 
the  output  from  that  run  was  compared  to  the  final  condition  given  in  the  original  data.  To 
facilitate  this  comparison,  the  Surfer  software  was  used  to  estimate  both  the  initial  and  final 
water  table  elevation  and  contaminant  concentration  for  each  finite  difference  node.  The  RMS 
error  of  the  differences  between  the  heads  or  concentrations  over  the  entire  finite  difference 
grid  was  calculated  for  a  given  run  of  Bioplume-II  using  the  FORTRAN  program  listed  in 
Appendix  D. 

4.2.3  One-Factor-at-a-Time  Approach.  The  one-factor-at-a-time  calibration  con¬ 
sisted  of  two  phases:  the  flow  calibration  and  the  transport  calibration. 

4.2.3. 1  One-Factor-at-a-Time  Flow  Calibration.  The  first  parameter  varied, 
chosen  arbitrarily,  was  the  transmissivity.  Starting  with  the  nominal  value  of  0.00226042 
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mentioned  in  paragraph  4.2. 1.4,  the  transmissivity  was  incrementally  varied.  The  goal  was  to 
minimize  the  error  criterion  by  varying  the  value,  within  a  reasonable  range. 

Table  1  shows  the  values  used  and  the  resulting  error  values.  The  range  of  transmis¬ 
sivities  considered  was  based  on  the  feasible  range  of  hydraulic  conductivities  in  the  United 
States  (56:8).  Walton  lists  0.001  ft  per  day  through  10,000  ft  per  day  as  the  range  for  hy¬ 
draulic  conductivity  in  the  United  States.  The  range  for  the  transmissivity  values  was  based 
on  the  fact  that  transmissivity  can  be  expressed  as  the  product  of  aquifer  depth  and  hydraulic 
conductivity,  as  depicted  in  equation  14.  The  aquifer  was  modeled  as  1.97  feet  deep,  which 
was  the  depth  of  the  screened  portion  of  the  injection  well.  Therefore  the  transmissivity  was 
1.97  times  the  hydraulic  conductivity.  Transmissivity  variation  was  bounded  between  2E-8 
and  0.228  ft^/sec,  as  depicted  in  equations  14  and  15. 


Run  No. 

Transmissivity  (ft/sec) 

RMSflow 

1 

2.26042E-3 

2.23746 

2 

3.26042E-3 

2.23746 

3 

1.26042E-3 

2.23746 

4 

2.26042E-2 

2.23746 

5 

2.26042E-4 

2.23746 

6 

2.26042E-1 

2.23746 

7 

2.26042E-5 

2.23746 

8 

2.26042E-6 

2.23746 

9 

2.26042E-7 

2.23746 

10 

2.26042E-8 

2.23746 

Table  1.  Values  of  Transmissivity  and  RMS  Error 


T  =  Kb 


(12) 


/0.001ft\ 
V  day  ) 


f  Iday  \ 

V86400sec; 


(1.97ft)  =  2*  10-®ftVsec 


(13) 
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10000ft 

day 


(14) 


I  f  1  (1.97ft)  =  0.228ftVsec 

J  V86400sec  /  ^  ^ 

Table  1  shows  that  increasing  and  decreasing  the  transmissivity  by  l.OE-3  from  its 
initial  value  of  2.26042E-3  did  not  alter  the  error  criterion  for  the  flow  calibration.  The  value 
remained  at  2.23746  feet,  regardless  of  the  value  of  the  transmissivity  parameter.  Scaling 
the  transmissivity  up  and  down  by  powers  of  ten  also  did  not  change  the  value  of  the  error 
criterion  for  the  flow  calibration. 

The  second  parameter  varied  was  porosity.  Boggs  and  others  describe  the  aquifer  as 
"sandy  gravel  and  gravelly  sand"  (10:3282).  Walton  lists  a  range  of  reasonable  permeabilities 
for  a  sand  and  gravel  aquifer  as  0.20  through  0.35  (56:413).  Boggs  and  others  (10:3282)  took 
samples  at  the  site  and  determined  that  the  porosities  at  the  site  had  either  a  mean  of  0.30  with 
a  standard  deviation  of  0.07  or  a  mean  of  0.32  with  a  standard  deviation  of  0.09.  To  include 
a  wide  range  of  reasonable  values  the  porosity  calibration  was  bounded  between  0.20  and 
0.41.  The  results  of  the  porosity  calibration  are  shown  in  Table  2.  These  results  show  that  the 
variation  of  porosity  within  a  reasonable  range  did  not  reduce  the  error  criterion. 


Run  No. 

Porosity 

RMS  flow 

11 

0.215 

2.23746 

12 

0.410 

2.23746 

Table  2.  Values  of  Porosity  and  RMS  Error  Criteria 

The  third  parameter  varied  was  the  recharge  rate.  There  was  no  data  available  on  a 
reasonable  range  for  the  recharge  rate,  so  the  model  was  started  with  a  few  values  to  determine 
which  might  yield  a  good  initial  value.  Simulation  time  was  used  to  determine  a  feasible 
starting  value  for  this  parameter.  To  predict  how  much  CPU  time  the  simulation  would  require 
for  a  given  run,  the  number  of  particle  moves  required  per  time-step  was  recorded  for  several 
runs.  Any  value  more  than  approximately  50  moves  was  prohibitive  in  terms  of  time  to  run 
and  memory  requirements.  Therefore  the  recharge  rate  was  changed  from  zero  to  5E-7  feet 
per  second.  Refer  to  the  Bioplume  II  manual  for  more  information  on  particles  in  Bioplume 
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II  (46:7-6).  Table  3  shows  the  number  of  particle  moves  per  time  step  required  by  Bioplume 
II  for  several  diffuse  recharge  rates  as  displayed  by  Bioplume  II. 


Diffuse  Recharge  in  feet  per  second 

Required  Particle  Moves  Per  Time  Step 

Estimated  Run  Time 

5E-3 

58086 

5E-4 

5809 

5E-5 

581 

3  days 

5E-7 

7 

1  hour 

Table  3.  Diffuse  Recharge  and  Required  Particle  Moves 


A  somewhat  arbitrary  range  of  -lE-6  through  5E-6  was  established,  based  on  computer 
limitations  described  above.  Results  of  the  calibration  using  diffuse  recharge  are  shown  in 
Table  4.  Note  that  recharge  is  negative  discharge.  The  results  are  also  displayed  in  graphical 
form  in  Figure  1. 


Run  No. 

Diffuse  Discharge  in  feet  per  second 

RMS  flow 

13 

-5E-7 

4.54487 

14 

5E-7 

0.990886 

15 

lE-6 

2.92961 

16 

7.5E-7 

1.82458 

17 

2.5E-7 

120544 

Table  4.  Diffuse  Discharge  and  RMS  Error  Criteria 


Clearly,  the  hydraulic  heads  are  most  sensitive  to  the  recharge/discharge  parameter  when 
compared  to  the  other  parameters  studied.  The  discharge  was  set  to  5E-7  feet  per  second, 
and  since  all  three  flow  parameters  had  already  been  tried,  transmissivity  was  randomly 
selected  to  be  recalibrated.  During  the  recalibration  of  transmissivity,  its  variation  affected 
response.  These  results  are  shown  on  Table  5  and  displayed  graphically  in  Figure  2.  The  best 
transmissivity  parameter  value  selected  was  2.76042E-3  feet/second  with  an  RMS  value  of 
0.921225  feet. 

Because  the  transmissivity  recalibration  did  not  bring  the  error  criterion  below  a  desired 
cutoff  value  of  0.5  feet  the  calibration  was  continued.  Porosity  was  randomly  selected  for 
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Possible  Singularity  at  Discharge=2.5E-7  ft/sec 

0.5 1 - 1 - 1 - 1 - ■ - 1 - 1 - 1 - 1 - 

-6  -4  -2  0  2  4  6  8  10 

Discharge  (ft/sec)  ^  ^  q-7 

Figure  1.  RMS  Error  vs.  Diffuse  Discharge  Rate 

recalibration.  The  results  of  the  porosity  recalibration  are  shown  in  Table  6.  Variation  of 
porosity  once  again  did  not  affect  the  error  criteria. 

Because  the  range  of  reasonable  porosity  values  does  not  permit  expansion  by  a  factor 
of  ten  in  either  direction,  recalibration  on  porosity  was  complete.  Discharge  was  randomly 
selected  for  recalibration  by  assigning  ranges  to  each  factor  and  using  a  random  number 
generator. 


Run  No. 

Transmissivity  (ft/sec) 

RMS  flow 

18 

3.26042E-3 

0.990044 

19 

4.26042E-3 

1.19303 

20 

3.76042E-3 

1.09244 

21 

2.76042E-3 

0.921225 

22 

2.51042E-3 

0.928127 

23 

3.01042E-3 

0.947168 

Table  5.  Recalibration  on  Transmissivity 


2.5  3  3.5  4  4.5 


T  ransmissivity  (ft''2/sec)  x  1 0'* 

Figure  2.  RMS  Error  vs.  Transmissivity 


Table  6.  Recalibration  on  Porosity 


The  discharge  recalibration  began  with  the  current  5E-7  feet  per  second  and  moved  a 
step  in  both  directions.  Because  no  improvement  was  identified  in  the  response  an  order  of 
magnitude  jump  would  have  been  taken  in  both  directions.  Due  to  CPU  time  and  memory 
restrictions  the  jump  could  only  be  taken  in  the  direction  of  smaller  values.  The  largest  value 
used  was  lE-6  feet  per  second. 

The  results  of  the  discharge  recalibration  are  shown  in  Table  7  and  displayed  graphically 
in  Figure  3.  Even  with  this  jump  no  improvement  was  found.  Due  to  lack  of  progress  the 
one-factor-at-a-time  flow  calibration  was  determined  to  be  completed  at  this  point  and  the 
transport  calibration  begun. 
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Run  No. 

Diffuse  Discharge  (ft/sec) 

RMS  flow 

25 

5E-7 

0.921225 

26 

lE-6 

2.11809 

27 

0.0 

2.23746 

28 

5E-8 

2.05905 

Table  7.  Recalibration  on  Discharge 
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2 
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1.4 
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With  the  completion  of  calibration  by  transmissivity,  porosity  and  diffuse  recharge,  the 
flow  calibration  was  completed.  In  a  total  of  28  runs  the  RMS  flow  criteria  was  reduced  to 
0.921225  feet. 

4.2.3.2  One-Factor-at-a-Time  Transport  Calibration.  After  the  flow  cali¬ 
bration  was  completed,  the  transport  calibration  was  begun.  The  first  parameter  varied  was 
the  longitudinal  dispersivity.  Hess  cites  Young  and  Boggs  (64:2012)  who  performed  field 
experiments  on  the  MADE-2  site  to  approximate  the  longitudinal  dispersivity  to  be  between 
10  and  42  meters  (32.8  and  137.76  feet).  Using  these  values  as  the  reasonable  range,  the 
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Discharge  (ft/sec)  ^  0-7 

Figure  3.  RMS  Error  vs.  Discharge 
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calibration  was  continued.  Table  8  presents  the  results,  which  show  that  the  variation  of 
longitudinal  dispersivity  did  not  change  the  error  criteria. 


Run  No. 

Longitudinal  Dispersivity 

RMS  transport 

1 

46.3 

67.1831 

2 

32.8 

67.1831 

3 

56.3 

67.1831 

4 

137.76 

67.1831 

Table  8.  Values  of  Longitudinal  Dispersivity  and  RMS  Error 


The  second  parameter  varied  was  the  ratio  of  longitudinal  to  transverse  dispersivity. 
This  value  had  been  held  at  0.4,  but  was  then  varied  between  the  limits  suggested  in  the 
Bioplume  manual  (46:7-15)  of  0.001  through  1.0.  The  results  of  calibration  on  this  ratio  are 
shown  in  Table  9  and  demonstrate  that  varying  the  value  of  transverse  dispersivity  did  not 
change  the  response  of  the  model. 


Run  No. 

Ratio  of  Longitudinal  to  Transverse  Dispersivity 

RMS  transport 

5 

0.5 

67.1831 

6 

0.3 

67.1831 

1 

1.0 

67.1831 

8 

0.04 

67.1831 

9 

0.004 

67.1831 

Table  9.  Ratio  of  Longitudinal  to  Transverse  Dispersivity  and  RMS  Error 


The  third  parameter  varied  was  the  stoichiometric  ratio  for  reaction  between  oxygen  and 
hydrocarbon  under  aerobic  biodegradation.  This  value  is  the  ratio  of  oxygen  to  carbon  atoms 
required  for  biodegradation  of  benzene.  In  the  case  of  complete  mineralization  of  benzene  to 
carbon  dioxide  and  water,  the  ratio  may  be  calculated  exactly  from  the  stoichiometric  equation 
for  aerobic  mineralization  of  benzene  shown  below  (48:56).  Note  that  fifteen  oxygen  atoms 
are  needed  for  every  six  carbon  atoms.  The  stoichiometric  ratio  is  therefore  15/6  or  2.5.  The 
value  must  be  calibrated  because  full  mineralization  of  the  benzene  may  not  occur  over  the 
modeled  time  period.  This  value  was  varied  between  2  and  4  based  on  the  recommendation 
of  Dr.  Rifai  (45).  The  results  of  calibration  on  the  stoichiometric  ratio  are  shown  in  Table  10. 


Run  No. 

Stoichiometric  Ratio 

RMS  transport 

10 

4.0 

67.1831 

11 

2.0 

67.1831 

Table  10.  Stoichiometric  Ratio  and  RMS  Error 


The  fourth  parameter  to  be  varied  was  the  Retardation  Factor.  This  factor  represents 
the  delay  in  spreading  of  the  contaminant  plume  caused  by  adhesion  to  the  soil  matrix  and 
is  defined  as  shown  in  equation  17.  The  value  of  0w,  percent  saturation,  is  unity  in  this 
case  because  the  simulation  models  the  saturated  zone.  Simplification  results  in  equation  19. 
Before  this  point,  the  Retardation  Factor  was  given  a  value  of  one.  This  assumption,  in  effect, 
stated  that  Kd,  the  distribution  coefficient  of  the  aquifer  or  pB,  the  bulk  density  of  the  solid 
matrix,  was  zero. 


Rf  —  1  + 


PbKd 

6w 


(15) 


i?/  =  1  +  pbKd  (16) 

Kd  =  (17) 

PB 

Stauffer  and  others  (50:75)  list  a  mean  value  of  retardation  of  1.20  for  the  site  with 
a  standard  deviation  of  0.20.  Arbitrarily  taking  plus  or  minus  one  standard  deviation  as 
a  reasonable  range  for  the  factor,  retardation  was  allowed  to  vary  between  1.00  and  1.40. 
Because  Bioplume  treats  retardation  as  a  function  of  pB  and  Kd  this  range  had  to  be  reflected 
in  their  values.  Bulk  density  was  set  to  1.77  grams  per  cubic  centimeter,  a  mean  value 
from  field  data  (50:5).  Equation  19  was  then  solved  for  Kd  as  shown  in  equation  20.  This 
equation  was  then  used  to  find  a  range  for  the  Distribution  Coefficient  given  the  range  for  the 
Retardation  Factor.  The  range  selected  for  the  Distribution  Coefficient  was  [0.00, 0.22599] 
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with  units  of  cubic  centimeters  per  gram.  The  results  of  the  variation  of  retardation  factor 
were  uninteresting.  Five  values  of  Kd  were  tried  with  the  same  result  for  each  run  as  shown 
in  Table  11. 


Run  No. 

Distribution  Coefficient 

RMS  transport 

12 

0.11299 

67.1831 

13 

0.21299 

67.1831 

14 

0.01299 

67.1831 

15 

0.22599 

67.1831 

16 

0.01299 

67.1831 

Table  11.  Distribution  Coefficient  and  RMS  Error 


The  last  transport  calibration  factor  to  be  tried  was  the  reaeration  coefficient.  Five 
values  of  this  factor  were  used  as  shown  in  Table  12. 


Run  No. 

Reaeration  Coefficient 

RMS  transport 

17 

1.184E-3 

67.1831 

18 

2.184E-3 

67.1831 

19 

0.184E-3 

67.1831 

20 

5E-3 

67.1831 

21 

1.184E-4 

67.1831 

Table  12.  Reaeration  coefficient  and  RMS  Error 


Once  again,  variation  of  this  factor  did  not  lead  to  variation  in  the  response.  Each  time 
the  model  was  run,  the  results  were  the  same:  RMS  transport  was  67.1831.  Due  to  lack  of 
progress  the  one-factor-at-a-time  flow  calibration  was  determined  completed,  having  made  no 
progress  after  21  runs  of  the  model. 

4.2.4  RSM  Approach.  The  second  calibration  procedure  used  was  the  Response 
Surface  Methodology  technique.  This  procedure  involved  three  distinct  phases.  The  screening 
phase  is  designed  to  determine  if  any  factors  are  insignificant  and  may  be  eliminated  from 
further  consideration.  The  first-order  phase  uses  an  empirical  first-order  regression  model  to 
approximate  the  true  response  surface  and  find  parameter  combinations  that  lead  to  optimal 
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response  values.  The  final  phase  is  the  second-order  phase  in  which  second-order  regression 
models  are  used  to  improve  the  empirical  model  and  find  optimal  parameter  combinations. 

4.2.4. 1  RSM  Flow  Calibration.  The  RSM  flow  calibration  consisted  of  the 
three  phases  described  above. 

Screening  Phase  of  RSM  Flow  Calibration.  The  screening  design 
included  three  factors:  the  transmissivity,  porosity  and  diffuse  discharge.  A  full  factorial 
design  was  used  because  of  its  simplicity.  Eight  runs  (2^)  runs  of  Bioplume  were  used  to 
determine  which  factors,  if  any,  could  be  excluded  from  the  analysis. 

These  runs  were  performed  at  all  possible  combinations  of  high  and  low  levels  of  the 
range  of  each  parameter.  The  levels  are  shown  in  Table  13.  Due  to  the  cost  in  terms  of  CPU 
time  and  memory,  the  discharge  could  not  be  set  with  an  absolute  value  higher  than  lE-6. 
Transmissivity  settings  above  0.0228  also  caused  problems  with  CPU  time  and  memory.  These 
values  were  high  enough,  however,  to  provide  a  sufficient  range  for  the  screening  design.  The 
results  of  the  eight  runs  are  shown  in  Table  14. 


Factor 

Low-Level 

High-Level 

Transmissivity 

2E-8 

0.0228 

Porosity 

0.20 

0.41 

Discharge 

-lE-6 

lE-6 

Table  13.  Levels  of  Flow  Parameters 


Once  the  runs  were  completed  they  were  evaluated  using  regression.  A  stepwise 
procedure  adds  and  deletes  one  factor  or  interaction  term  at  a  time  in  the  regression  equation 
to  determine  whether  the  term  improves  the  fit  of  the  regression.  Only  those  terms  significantly 
improving  the  fit  the  are  included  at  each  step.  This  procedure  was  used  to  determine  which 
factors,  including  interactions,  were  worth  keeping  in  the  analysis.  Since  the  Bioplume-II 
model  produces  deterministic  output,  the  statistical  significance  of  the  regression  model  or  of 
any  individual  factors  was  not  meaningful.  SAS  was,  however,  still  able  to  provide  a  least 
squares  model  of  the  error  response.  The  stepwise  procedure  approximated  a  parsimonious 
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Run  No. 

Transmissivity 

Porosity 

Discharge 

RMS  flow 

1 

2E-8 

0.20 

-lE-6 

27246.5 

2 

2E-8 

0.20 

lE-6 

544968 

3 

2E-8 

0.41 

-lE-6 

27246.5 

4 

2E-8 

0.41 

lE-6 

544968 

5 

0.0228 

0.20 

-lE-6 

2.22 

6 

0.0228 

0.20 

lE-6 

1.81 

7 

0.0228 

0.41 

-lE-6 

2.22 

8 

0.0228 

0.41 

lE-6 

1.81 

Table  14.  Screening  Design  and  RMS  Error 


model,  using  the  coded,  i.e.  transformed,  data  rather  than  the  actual  parameter  values.  The 
resulting  regression  model  revealed  that  neither  porosity  nor  any  of  the  interaction  terms  that 
included  porosity  were  worth  including  in  the  calibration.  The  stepwise  procedure  from  SAS 
is  shown  in  Appendix  C-1. 


First-Order  Design  Phase  ofRSM  Flow  Calibration.  The  First-Order 
design  phase  ofRSM  was  used  with  only  the  two  significant  factors  from  the  screening  design: 
transmissivity  and  discharge.  Given  that  a  full  factorial,  2^  design  was  used  for  these  two 
factors,  four  runs  of  Bioplume  were  needed  for  each  design  used.  The  high  and  low  settings 
for  the  first-order  design  phase  differed  from  those  of  the  screening  phase.  In  the  screening 
phase,  the  settings  spanned  the  range  of  reasonable  values  for  the  variable.  In  the  first-order 
design  phase,  they  represented  a  range  around  the  current  setting.  A  range  of  ±80%  of  the 
initial  value  was  used  (19:4-13).  The  results  of  the  first  run  of  the  first-order  design  phase  are 
shown  in  Table  15. 


Run  No. 

Transmissivity 

Discharge 

RMS  flow 

9 

4.068756E-3 

5E-8 

2.11605 

10 

4.068756E-3 

-5E-8 

2.36023 

11 

4.52084E-4 

5E-8 

1.24133 

12 

4.52084E-4 

-5E-8 

3.37311 

Table  15.  First  Run  of  the  First-Order  Design  Phase  and  RMS  Error 
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These  results  were  modeled  in  a  first-order  linear  regression  equation  using  SAS.  The 
regression  output  from  SAS  is  shown  in  Appendix  C-2.  The  SAS  output  shows  that  the  fit 
was  poor,  with  an  value  of  0.6139.  In  order  to  develop  a  more  reliable  model  (that  could 
be  trusted)  the  design  size  was  reduced  to  ±40%  of  the  factor’s  value.  A  centerpoint  value  of 
8.401E-7  was  randomly  selected  for  discharge  because  a  range  around  a  value  of  zero  could 
not  be  based  on  a  percentage.  The  results  of  runs  from  this  design  are  shown  in  Table  16. 


Run  No. 

Transmissivity 

Discharge 

RMS  flow 

13 

1.356252E-3 

5.0406E-7 

i.imi 

14 

1.356252E-3 

1.176E-6 

7.46751 

15 

3.164588E-3 

5.0406E-7 

0.967852 

16 

3.164588E-3 

1.176E-6 

2.21034 

Table  16.  Second  Run  of  the  First-Order  Design  Phase  and  RMS  Error 


These  results  were  also  fit  to  a  first-order  regression  equation  using  SAS.  The  value 
was  an  improved  0.8398,  but  the  F-ratio  was  2.622.  The  F-test  has  no  meaning  here  because 
it  is  based  on  an  assumed  normal  error  distribution  around  each  response  value.  Although  the 
response  surface  is  an  error  surface,  the  error  is  deterministic  and  so  the  normal  assumption 
can  not  hold.  The  F-ratio  is  the  amount  of  the  total  sum  of  squared  response  which  is  explained 
by  the  model  divided  by  that  attributed  to  lack  of  fit.  Therefore  any  value  over  1 .0  means  more 
of  the  sum  of  squares  is  explained  by  the  model  than  is  attributed  to  lack  of  fit.  The  higher  the 
F-ratio,  the  better  the  fit.  The  regression  output  is  shown  in  Appendix  C-3. 

The  second  attempt  to  develop  a  significant  model  was  the  addition  of  a  center  point 
to  the  2^  design.  This  made  for  a  total  of  five  data  points.  The  fifth  data  point  is  shown  in 
Table  17. 


Run  No. 

Transmissivity 

Discharge 

RMS  flow 

17 

2.26042E-3 

8.401E-7 

2.21011 

Table  17.  Centerpoint  of  Previous  Design  and  RMS  Error 
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The  fifth  data  point  did  not  improve  the  value,  but  it  did  increase  the  F-ratio  to  4.370. 
This  was  thought  to  be  acceptable  for  developing  a  gradient  vector.  The  regression  output  is 
shown  in  Appendix  C-4.  The  coefficients  of  the  regression  equation  were  used  to  calculate  the 
elements  of  the  gradient  vector.  The  vector  desired  was  the  vector  of  steepest  descent,  which 
is  based  on  the  additive  inverses  of  each  element.  The  elements  had  to  first  be  normalized  to 
establish  the  unit  gradient  vector.  Then  the  values  were  decoded  because  they  were  based  on 
the  transformed  coordinate  system,  where  the  center  of  the  current  design  had  coordinates  of 
(0.0). 

The  transformation  followed  equation  4,  which  was  solved  for  9k,  as  shown  in  equa¬ 
tion  20.  Note  that  is  the  coefficient  value,  in  coded  units,  from  the  regression  equation. 
Sk  is  the  half  range  of  the  region  in  standard  units  and  6ko  is  the  value  of  the  parameter,  in 
standard  units,  at  the  center  of  the  design  region.  The  results  are  shown  in  Table  18. 


9k  =  XkSk  +  9ko 


(18) 


Factor 

Regression  Coefficient 

Normalized  Coefficient 

Uncoded  Value 

Transmissivity 

-1.625022 

-0.707154 

1.6210E-3 

Discharge 

1.624807 

0.707060 

1.0777E-6 

Table  18.  Table  of  Calculations  for  Steepest  Descent  Gradient 


The  uncoded  values,  which  were  the  elements  of  the  uncoded  steepest  ascent  gradient 
vector,  were  then  subtracted  from  the  parameter  values  at  the  center  of  the  design  region. 
The  new  coordinate  represented  a  single  step  along  the  gradient  vector.  A  new  test  at  this 
point  revealed  that  the  design  region  had  been  too  large.  The  response  at  the  unit  step  was 
worse  than  before,  when  it  was  expected  to  to  have  generated  an  improved  response.  The  new 
response  is  shown  in  Table  19  as  run  number  18. 

Using  the  same  design,  the  step  size  was  decreased  in  the  direction  of  the  gradient 
vector.  Instead  of  moving  one  unit  vector  in  the  gradient  direction,  a  new  test  point  was 
established  at  50%  of  the  steepest  descent  gradient  vector.  This  test  point  revealed  a  response 
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of  0.941898  for  flow  calibration.  This  value  was  an  improvement  over  the  response  at  the 
center  of  the  first  design,  2.21077  feet.  Values  at  40  and  60%  of  the  unit  step  were  taken  to 
determine  whether  the  50%  step  led  to  the  best  response  along  the  gradient  vector. 


Run  No. 

Step  Number 

RMS  Flow 

17 

Center  of  First  Design 

i.imi 

19 

4/10 

1.16225 

20 

5/10 

0.941898 

21 

6/10 

1.01111 

18 

Unit  Step 

6.15908 

Table  19.  Table  of  RMS  Flow  Values  along  First  Gradient  Vector 


4.5  5  5,5 

Step  Number 


6.5 


Figure  4.  Graph  of  RMS  Flow  Along  Gradient  Vector 


Once  the  50%  step  was  chosen  as  the  low  point  along  the  descent  path  a  new  design 
was  established.  The  first  design,  centered  on  the  50%  step,  that  led  to  an  F-ratio  over  two 
was  at  ±5%  of  the  parameter  values.  These  results  are  shown  in  Table  20.  It  had  an  F-ratio 
of  5.821  and  an  value  of  85.34%,  see  Appendix  C-5  for  the  SAS  output  from  this  design. 
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Run  No. 

Transmissivity 

Discharge 

RMS  flow 

22 

1.39E-3 

2.86E-7 

0.940928 

23 

1.39E-3 

3.16E-7 

1.01783 

24 

1.53E-3 

2.86E-7 

0.918875 

25 

1.53E-3 

3.16E-7 

0.942795 

Table  20.  Run  of  a  5%  Design  Centered  on  5th  Step 

The  next  step  was  to  calculate  the  steepest  descent  gradient  from  this  point.  Using 
equation  20,  the  steepest  descent  gradient  was  calculated  from  the  coefficients  of  the  regression 
equation.  The  results  of  these  calculations  are  shown  in  Table  21 . 

The  unit  gradient  step  was  not  the  first  step  used  for  the  gradient  search  here  because 
the  unit  gradient  step  would  have  been  beyond  the  feasible  region  for  transmissivity.  Instead, 
in  an  attempt  to  explore  the  first-order  design  further,  10%  steps  were  taken.  Then  step  ten 
was  taken  in  the  direction  of  the  discharge  component  of  the  steepest  descent  gradient  vector. 
This  tenth  step  is  the  gradient  search  subject  to  a  constraint  as  outlined  in  Box  and  Draper 
(12).  The  results  of  the  Bioplume  II  runs  are  shown  in  Table  22.  Note  that  the  unit  gradient 
step  was  taken  at  step  nine  was  the  last  step  that  could  fit  within  the  range  of  acceptable  values 
for  transmissivity.  See  Figure  5  for  a  graph  of  the  progress  made  along  the  gradient  vector. 


Factor 

Regression  Coefficient 

Normalized  Coefficient 

Uncoded  Value 

Transmissivity 

-0.024272 

0.6936 

1.51  lE-3 

Discharge 

0.025206 

0.7203 

3.118E-7 

Table  21.  Table  of  Calculations  for  Second  Steepest  Descent  Gradient 

Finally  a  5%  design  around  the  ninth  step  along  the  last  gradient  vector  was  tried.  This 
design  had  results  shown  in  Table  23.  The  SAS  analysis  of  this  first-order  design  showed  an 
F-ratio  of  1.840  and  the  value  was  only  0.6479.  See  Appendix  C-6  for  these  results. 

Second-Order  Design  Phase  ofRSM  Flow  Calibration.  Because  this 
design  was  already  so  small,  a  second-order  design  was  appropriate  at  this  stage.  Again  the 
reader  should  consult  Box  and  Draper  (12)  for  more  information  on  second-order  strategies. 
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Table  22.  Table  of  RMS  How  Values  along  Second  Gradient  Vector 


Run  No. 

Transmissivity 

Discharge 

RMS  flow 

36 

9.823E-5 

1.9361E-8 

0.925219 

37 

9.823E-5 

2.1399E-8 

0.977799 

38 

1.0857E-4 

1.9361E-8 

0.923916 

39 

1.0857E-4 

2.1399E-8 

0.925219 

Table  23.  Run  of  a  5%  Design  Centered  on  9th  Step 


Four  axial  points  were  added  at  ±a  along  each  of  the  two  axes  to  develop  a  central  composite 
design,  a  was  chosen  to  equal  the  fourth  root  of  the  number  of  factorial  points  in  the  full 
factorial  design,  1.414  in  this  case,  as  suggested  by  Auclair  (5).  The  additional  axial  points 
are  shown  in  Table  24. 


Run  No. 

Transmissivity 

Discharge 

RMS  flow 

40 

1.1071E-4 

2.038E-8 

0.919488 

41 

9.60885E-5 

2.038E-8 

0.958185 

42 

1.034E-4 

2.1821E-8 

0.954969 

43 

1.034E-4 

1.8939E-8 

0.919911 

Table  24.  Additional  Points  Centered  on  9th  Step 
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Figure  5.  Graph  of  RMS  Flow  Along  Second  Gradient  Vector 


All  nine  points  in  the  design  on  the  9th  step  were  then  run  through  a  second  order  model 
in  SAS.  This  model  approximated  the  coefficients  of  equation  21;  the  SAS  results  are  shown 
in  Appendix  C-7  and  include  a  canonical  and  a  ridge  analysis. 


y  —  ^0  +  ^tT  +  ^dD  +  ^ttT^  +  /3ddD^  +  I^tdTD  (19) 

The  canonical  analysis  shows  a  local  minimum  could  exist  at  a  transmissivity  value  of 
-0.595968  (1.00319E-4  uncoded)  and  a  discharge  value  of  -1.689911  (1.8658E-8  uncoded). 
This  point  was  run  through  Bioplume  for  an  RMS  flow  value  of  0.918963  feet. 

The  ridge  analysis  used  SAS  to  extrapolate  out  from  the  design  region  to  find  points 
that  would  be  good  candidates  for  further  exploration.  The  results  of  the  ridge  analysis  are 
shown  in  Appendix  C-7.  Four  points  along  the  predicted  ridge  were  tested  to  see  if  a  ridge  of 
decreasing  response  really  existed  and  if  the  responses  were  better  than  at  the  stationary  point 
predicted  by  the  canonical  analysis.  See  Table  25  for  results  of  the  ridge  analysis.  Note  that 
the  distance  along  the  ridge  given  in  the  table  is  in  coded  units. 
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Run  No. 

Distance  along  Ridge 

Transmissivity  (coded) 

Discharge  (coded) 

RMS  flow 

44 

1 

1.03696E-4  (0.057195) 

2.13973E-8  (-0.998363) 

0.942336 

45 

2 

9.95182E-5  (-0.750829) 

1.84911  (-1.853714) 

0.918993 

46 

3 

9.5795E-5  (-1.470982) 

1.77157  (-2.614615) 

0.919196 

A1 

4 

9.2159E-5  (-2.174265) 

1.69587  (-3.357465) 

0.919517 

Table  25.  Ridge  Analysis  for  RSM  Flow  Calibration 


Because  the  ridge  analysis  did  not  improve  the  response  beyond  the  value  of  the 
stationary  point  the  calibration  was  stopped.  The  best  response  found  during  the  ridge  analysis 
was  0.918993  feet  at  a  transmissivity  of  9.95182E-5  square  feet/second  and  a  discharge  of 
1 .849 1 1  feet/second. 

The  best  response  found  was  not  actually  at  the  stationary  point,  but  at  the  corner  of  the 
second  of  the  first  order  designs.  This  value  occurred  at  a  transmissivity  of  1.53E-3  square 
feet  per  second  and  a  discharge  of  2.86E-7  feet  per  second.  The  root  mean  squared  error 
criterion  at  this  point  was  0.918875  feet,  which  was  less  than  a  ten  thousandth  of  a  foot  lower 
in  response  than  the  stationary  point.  If  RSM  had  been  halted  at  this  point,  then  a  superior 
response  would  have  been  obtained  using  RSM  in  only  24  runs  of  the  Bioplume-II  model. 
The  RSM  calibration  was  not  terminated,  however,  because  in  exploring  the  capabilities  of 
the  procedure  the  specific  response  at  run  number  24  was  simply  overlooked. 

4.2.4.2  RSM  Transport  Calibration.  For  the  transport  calibration  five  factors 
had  to  be  considered.  These  were  longitudinal  dispersivity  {ocl),  the  ratio  of  transverse  to 
longitudinal  dispersivity  (^),the  stoichiometric  ratio  (F),  the  retardation  factor  {Rj),  and  the 
Reaeration  Coefficient  (RC). 

Screening  Phase  of  RSM  Transport  Calibration.  Using  a  full  factorial 
design,  a  screening  phase  was  run  for  all  the  factors.  This  design  included  2^  =  32  runs  at  the 
high  and  a  low  levels  for  each  factor.  These  levels  are  shown  in  Table  26,  and  the  results  of 
these  runs  are  given  in  the  SAS  output  in  Appendix  C-8. 
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Factor 

Low 

High 

OIL 

32.8 

137.46 

QTT 

Oil. 

0.001 

1.0 

F  (Stoichiometric  Coefficient) 

2.0 

4.0 

Kd  (Rf) 

0(1.00) 

0.11299(1.40) 

Reaeration  Coefficient 

0 

0.005 

Table  26.  Details  of  Screening  Design  Levels 


Although  the  screening  design  suggested  that  the  stoichiometric  ratio  and  the  retardation 
factor  had  the  least  effect  on  the  response,  the  response  surface  required  special  treatment. 
Two  of  the  32  runs  had  response  values  below  the  67.1831  that  resulted  from  most  of  the 
One-Factor  calibration  approach.  These  values  were  obtained  only  when  the  stoichiometric 
ratio  was  set  to  a  high  value  and  the  Reaeration  coefficient  was  set  to  zero  (uncoded). 

Only  two  factors  were  considered  for  further  calibration,  the  longitudinal  dispersivity 
and  the  ratio  of  transverse  to  longitudinal  dispersivity.  The  stoichiometric  ratio  and  the 
reaeration  coefficient  were  set  to  their  high  and  low  values  respectively.  Because  of  it’s 
apparently  minimal  impact  on  response,  the  retardation  factor  was  set  tp  the  one,  i.e.  no 
retardation. 


First-  Order  Phase  ofRSM  Transport  Calibration.  The  starting  location 
selected  for  the  first-order  design  coincided  with  the  lowest  response  observed  in  the  screening 
phase.  This  point  could  not  be  used  as  the  center  of  a  new  design,  however,  because  it  was 
on  the  boundary  of  reasonable  values  of  the  parameters.  Therefore,  an  arbitrary  design  was 
chosen  that  included  the  minimum  point  as  a  corner.  This  choice  meant  setting  the  stochastic 
ratio  to  4.0  (+1  coded)  and  the  reaeration  coefficient  to  zero.  The  retardation  factor  was  set 
to  zero  (-1  coded)  because  this  value  led  to  the  lowest  response  in  the  screening  design.  The 
remaining  two  factors  were  then  varied  and  the  results  were  as  shown  in  Table  27. 

These  responses  were  then  fit  to  a  first-order  regression  model  using  SAS.  The  output 
is  shown  in  Appendix  C-9.  These  results  show  that  a  very  poor  fit  was  obtained. 
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run  No. 

Ot'p 

OtJ. 

RMS  Transport 

33 

123.714 

0.001 

61.1906 

34 

123.714 

0.002 

67.0921 

35 

137.46 

0.001 

67.0505 

36 

137.46 

0.002 

67.1115 

Table  27.  First-Order  Design  for  RSM  Transport  Calibration 

Second-Order  Phase  of  RSM  Transport  Calibration.  A  second-order 
design  was  next  used  to  improve  the  calibration.  Additional  tests  were  made  at  1.414  coded 
units  away  from  the  design  center  along  each  axis.  These  coordinates  would  have  rendered 
some  tests  out  of  the  feasible  range  of  factor  values.  To  avoid  this  problem,  the  readings  were 
taken  at  one  coded  unit  along  the  axis  when  necessary.  A  face  centered  central  composite 
design  would  have  been  a  more  conventional  choice  for  this  design,  but  this  unusual  design 
worked  acceptably  well.  The  center  of  the  design  was  also  included  for  a  total  of  five  additional 
points.  The  results  are  shown  in  Table  28. 


Run  No. 

ai  (coded) 

(coded) 

RMS  Transport 

37 

0.0015  (0) 

67.2464 

38 

137.46  (+1) 

0.0015  (0) 

67.0764 

39 

130.587  (0) 

67.0327 

40 

130.587  (0) 

0.001  (-1) 

67.0394 

41 

130.587  (0) 

0.0015  (0) 

67.0362 

Table  28.  Second-Order  Design  for  RSM  Transport  Calibration 


A  ridge  analysis  was  performed  using  the  nine  data  points  tested.  The  SAS  output 
is  shown  in  Appendix  C-10.  The  analysis  predicted  a  ridge  pointing  in  the  direction  of 
decreasing  longitudinal  dispersivity  and  increasing  transverse  dispersivity.  First  the  stationary 
point  identified  by  the  canonical  analysis  was  tested  and  then  points  were  tested  along  the 
predicted  location  of  the  ridge.  As  shown  in  Table  29,  no  change  in  response  was  identified. 
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Run  No. 

units  along  ridge 

ctL  (coded) 

RMS  Transport 

42 

Stationary  Point 

131.66(0.155432) 

0.00161  (0.220709) 

67.1831 

43 

1 

129.92  (-0.096925) 

0.0021  (1.200923) 

67.1831 

44 

2 

128.27  (-0.336868) 

0.0026  (2.202779) 

67.1831 

45 

3 

126.65  (-0.572427) 

0.0031  (3.184661) 

67.1831 

46 

4 

125.041  (-0.806931) 

0.0036(4.161754) 

67.1831 

Al 

5 

123.43  (-1.041017) 

0.0041  (5.136951) 

67.1831 

Table  29.  Ridge  Analysis  for  RSM  Transport  Calibration 


At  this  point  the  calibration  was  stopped  due  to  the  lack  of  sufficient  improvement  in 
the  response.  The  RSM  transport  phase  took  47  runs  and  had  a  best  response  of  67.0327  parts 
per  billion,  identified  at  ai,  =  130.587  and  ^  =  0.002207. 


4.3  Summary 

Table  30  summarizes  the  results  of  the  two  calibrations  to  allow  for  a  comparison.  It 
can  be  seen  that  RSM  took  longer  to  converge,  but  always  led  to  better  results. 


Method 

Best  Result 

Total  Runs 

Expected  Runs 

Flow  Calibration 

One-Factor 

0.921225 

36 

RSM 

0.918875 

Al 

24-43 

Transport  Calibration 

One-Factor 

67.1831 

21 

RSM 

67.0327 

Al 

23 

Table  30.  Summary  of  Calibrations 


Note  the  "Expected  Runs"  column  lists  the  number  of  runs  that  would  have  been  required 
if  efficient  Plackett-Burman  or  Fractional  Factorial  designs  had  been  used  for  the  Screening 
design.  The  range  shown  for  RSM  flow  calibration  (24-43  runs)  reflects  the  fact  that  the  best 
value  for  this  calibration  was  obtained  on  the  24th  run.  Since  all  runs  after  that  were  not  as 
good,  the  calibration  could  have  been  reasonably  stopped  at  any  time  after  that  run. 
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The  fact  that  RSM  took  more  runs  than  the  one-factor-at-a-time  approach  in  the  flow 
calibration  may  be  partially  because  RSM  was  used  in  an  exploratory  manner.  If  the  RSM  flow 
calibration  had  been  halted  when  it  surpassed  the  one-factor-at-a-time  approach,  only  19  runs 
would  have  been  required  with  efficient  Plackett-Burman  or  Fractional  Factorial  Screening 
designs. 

The  comparison  for  the  transport  calibration  does  not  present  the  full  picture  either.  The 
one-factor-at-a-time  approach  made  no  progress  at  all.  To  take  the  time  to  find  the  a  response 
of  67.0327  feet  using  the  one-factor-at-a-time  approach  would  most  likely  have  taken  many 
more  runs  than  the  RSM  technique  took.  By  finding  an  area  of  improved  response  on  the 
transport  surface,  RSM  performed  better  than  the  one-factor-at-a-time  approach. 


\ 


44 


V.  Conclusions  and  Recommendations 


5.1  Conclusions 

The  results  were  not  impressive  for  either  technique.  This  was  likely  caused  by  the 
use  of  the  homogeneity  assumption  for  all  environmental  factors  calibrated.  The  MADE-2 
site  is  well  known  in  the  literature  for  its  heterogeneity  (50).  This  difference  between  the 
simulated  and  actual  conditions  could  prevent  significant  improvement  in  response  regardless 
of  calibration  technique. 

The  strength  of  RSM  could  have  become  more  apparent  if  zonation  had  been  used.  A 
zonated  model  would  allow  for  significant  improvement  in  response,  but  would  have  added 
significantly  to  the  number  of  factors  involved  in  the  calibration.  The  one-factor-at-a-time 
approach  would  have  been  tedious  to  use  for  a  zonated  model.  RSM  would  excel  here  due  to 
the  use  of  the  Screening  phase  and  the  application  of  efficient  Plackett-Burman  or  Fractional 
Factorial  designs.  RSM  would  save  runs  by  calibrating  only  the  most  significant  factors. 

The  major  weakness  of  the  one-factor-at-a-time  approach  demonstrated  by  this  research 
was  its  inability  to  take  interactions  between  parameters  into  consideration  during  the  calibra¬ 
tion  process.  By  varying  one  parameter  at  a  time  the  method  used  a  tunnel  vision  approach. 
In  effect,  the  approach  can  only  turn  at  right  angles  while  negotiating  a  response  surface. 
RSM  improves  on  this  approach  by  simultaneously  varying  every  parameter  to  improve  the 
response.  This  gradient  search  approach  allows  RSM  users  to  obtain  a  better  sense  of  the 
response  surface  and  move  in  any  direction  to  improve  the  response. 

One  difficulty  with  the  use  of  RSM  is  that  the  modeler  must  have  knowledge  of  some 
statistical  concepts  beyond  those  of  many  engineers  and  hydrogeologists.  A  solution  to  this 
problem  is  that  when  the  modeler  is  not  familiar  with  the  material,  another  individual  in  the 
organization  possibly  could  be  assigned  to  assist  with  the  calibration.  If  this  is  not  the  case, 
specialized  training  may  be  needed  to  enable  modelers  to  use  this  approach. 
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Another  problem  with  the  use  of  RSM  is  that  a  statistical  software  package  may  be 
needed  to  take  full  advantage  of  the  second-order  techniques.  Not  every  environmental 
engineering  organization  has  access  to  this  type  of  software.  However,  the  first-order  and 
some  of  the  second-order  techniques  may  be  used  with  a  spreadsheet  or  even  a  calculator. 
They  still  have  the  advantage  of  varying  more  than  one  factor  at  a  time  and  therefore  take 
interactions  between  parameters  into  consideration.  Therefore,  lack  of  access  to  an  advanced 
statistical  software  package  does  not  limit  the  modeler  to  the  one-factor-at-a-time  approach. 

RSM  took  longer  to  find  the  best  response  identified  than  the  one-factor-at-a-time 
approach.  This  was  partly  the  result  of  the  researcher’s  inexperience  with  the  method.  RSM 
also  took  longer  because  it  was  not  used  in  the  most  efficient  way  possible.  Full  factorial 
designs  were  used  all  the  way  through  the  research  because  of  their  simplicity.  If  Plackett- 
Burman  or  other  fractional  factorial  designs  had  been  used  the  number  of  runs  would  have 
been  reduced  significantly. 

RSM  showed  a  strength  in  the  calibration  of  the  transport  phase.  The  one-factor-at- 
a-time  approach  was  unable  to  make  any  progress  after  varying  all  five  transport  parameters 
throughout  a  feasible  range.  In  order  to  continue  the  calibration  using  this  approach,  a 
new  starting  point  would  be  selected  and  the  parameters  again  varied.  The  problem  is  that 
the  areas  of  improved  response  were  not  easy  to  find.  The  traditional  method  would  have 
difficulty  improving  upon  model  response  in  this  case.  Without  a  systematic  procedure  to 
vary  one  parameter  with  respect  to  the  others,  one-factor-at-a-time  calibration  is  little  more 
than  iterative  guessing. 

The  greatest  strength  of  RSM  is  in  its  ability  to  move  over  the  response  surface  while 
taking  interactions  between  parameters  into  account.  This  strength  allowed  it  to  find  a  region 
of  improved  response.  In  the  transport  case  the  improvement  was  roughly  0.1  ppb.  Although 
this  improvement  is  small,  it  is  real  because  the  model  is  deterministic.  The  region  was 
probably  a  narrow  and  shallow  depression  in  the  transport  calibration  error  response  surface, 
but  it  would  be  difficult  to  find  at  all  without  a  systematic  procedure.  RSM  found  it  easily  by 
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varying  all  the  parameters  simultaneously.  If  only  one  factor  is  varied  at  a  time,  many  features 
of  the  response  surface  will  not  be  identified. 

5.2  Recommendations 

For  further  research  in  this  area  the  following  are  recommended: 

1.  Reduce  the  simulated  time  from  15  months  to  a  few  months  in  order  to  avoid  the 
problem  created  by  the  plume  disappearing  in  the  simulation.  This  change  should 
reduce  the  prevalence  of  the  67. 1831  response  for  the  transport  calibration  and  improve 
the  behavior  of  the  transport  response  surface. 

2.  Zonate  some  of  the  parameters  to  increase  model  accuracy.  Zonation  of  parameters 
such  as  transmissivity  and  dispersivity  might  significantly  improve  model  response  and 
drop  the  average  error  response.  This  change  would  also  show  the  strength  of  RSM  in 
its  ability  to  vary  large  numbers  of  parameters  while  exploring  the  response  surface. 

3.  Repeat  the  exercise  on  another  data  set  or  using  another  model.  This  repetition  might 
provide  an  indication  of  how  RSM  compares  to  the  one-factor-at-a-time  approach 
overall,  rather  than  just  for  Bioplume  II  using  MADE-2  data,  as  in  this  research. 

4.  Calibrate  both  flow  and  transport  simultaneously  by  using  the  sum  of  the  RMS  flow  and 
transport  criteria  for  the  error  response  surface.  This  approach  could  improve  calibration 
by  taking  interactions  between  degree  of  flow  calibration  and  transport  calibration  into 
account. 

5.  Repeat  the  research  using  fractional  factorial  designs  to  see  how  far  the  number  of 
required  runs  could  be  reduced. 

6.  Repeat  the  research  using  several  different  error  criteria  to  compare  their  effect  on 
convergence. 

7.  Repeat  the  research  using  transient  hydraulics  to  see  whether  the  success  of  each 
calibration  technique  depends  on  that  aspect  of  the  modeling. 
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Appendix  A.  Biodegradation  of  Petroleum  Hydrocarbons 


The  contamination  of  soil  and  groundwater  has  become  a  major  issue  in  the  media  over 
the  past  20  years.  The  result  of  the  public’s  interest  in  the  issue  was  the  passage  of  laws  such 
as  the  Comprehensive  Environmental  Response,  Compensation  and  Liability  Act  (CERCLA) 
requiring  cleanup  of  contaminated  sites.  The  U.S.  Air  Force  necessarily  has  interest  in  the  issue 
due  to  the  number  of  contaminated  sites  located  on  Air  Force  (AF)  bases.  The  extensive  use  of 
jet  fuels  by  the  AF  has  led  to  a  special  interest  in  the  remediation  of  petroleum  hydrocarbons 
in  the  subsurface.  Due  to  the  high  solubility  of  some  aromatic  hydrocarbons,  such  as  benzene, 
toluene,  ethylbenzene,  and  xylenes  (BTEX),  they  pose  an  extra  risk  to  human  health  and  the 
environment.  This  appendix  will  therefore  focus  on  aromatic  hydrocarbons. 

Of  all  the  remediation  technologies,  natural  attenuation  stands  out  as  holding  unusual 
potential  to  decrease  the  cost  of  remediating  many  contaminated  sites.  Natural  attenuation  is 
a  remedial  action  where  subsurface  microorganisms,  in  combination  with  geophysical  factors, 
naturally  clean  the  soil  and  groundwater.  The  potential  savings  are  huge. 

The  natural  attenuation  process  involves  four  factors:  volatilization,  dispersion,  sorption 
and  degradation.  Volatilization  is  the  transfer  of  contamination  from  aqueous  to  vapor  phases. 
Dispersion  is  the  spreading  out  of  contanunants  that  reduces  concentration,  but  not  total 
mass.  Sorption  is  the  adsorption  of  contaminants  to  the  soil  matrix.  This  process  reduces 
concentration  of  contaminants  in  groundwater,  but  not  total  mass.  Degradation  includes  both 
abiotic  chemical  and  biotic  transformation  of  contaminant  species  into  other  substances.  This 
process  reduces  the  concentration  and  mass  of  contaminant.  McAllister  and  Chiang  (38:163- 
164)  notes  that  biodegradation  is  usually  the  most  significant  degradation  process  in  the 
subsurface.  Stauffer  and  others  (51:78)  supported  this  conclusion  by  showing  that  the  effect 
of  sorption  on  plume  stabilization  was  negligible  compared  to  the  effect  of  biodegradation  at 
a  site  under  aerobic  conditions. 

The  concept  of  biodegradation  of  petroleum  hydrocarbons  is  that  the  microorganisms 
use  the  contaminant  as  a  food  source.  The  general  equation  can  be  seen  as  (1 1:179): 
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Hydrocarbon  +  electron  acceptor  +  microorganisms  +  inorganic  nutrients 
Carbon  Dioxide  +  water  +  microorganisms  +  waste  products 

Biodegradation  is  classified  as  aerobic  if  the  electron  acceptor  is  molecular  oxygen 
and  anaerobic  otherwise.  Liss  and  Baker  (36:303)  list  nitrate,  sulfate,  carbonate,  iron  and 
manganese  oxide  as  alternative  electron  acceptors.  Aerobic  biodegradation  requires  at  least  1- 
2  ppm  dissolved  oxygen  and  can  rapidly  degrade  some  compounds.  First  order  rate  constants 
of  0.3- 1.3/ 

The  rate  of  biodegradation  may  be  measured  using  the  Respiration  Quotient  (30:477) 
equal  to  actual  respiration  over  potential  respiration,  where  potential  respiration  was  microbial 
oxygen  uptake  measured  after  addition  of  an  easily  biodegradable  carbon  source,  such  as  glu¬ 
cose.  They  found  a  correlation  coefficient  of  0.997  between  the  quotient  and  concentration  of 
polycyclic  aromatic  hydrocarbons  (PAH)  showing  the  respiration  coefficient  was  an  indicator 
of  microbial  metabolic  activity. 

Leahy  and  Colwell  (35:307)  note  the  impact  of  physical  factors  on  microbial  activity. 
They  cite  research  by  Atlas  and  Bartha  observing  that  low  temperatures  inhibit  microbial 
activity  and  research  by  Bossert  and  Bartha  observing  that  higher  temperatures  increase 
activity  up  to  a  point,  beyond  which  disruption  occurs.  They  also  cite  Dibble  and  Bartha 
in  concluding  that  extremely  high  concentrations  of  petroleum  hydrocarbons  may  decrease 
activity  due  to  toxicity  effects  of  components  in  the  petroleum  mixture. 

It  is  clear  that  some  toxicity  could  occur  because  of  the  wide  variety  of  compounds 
present  in  the  petroleum  hydrocarbon  mixture.  Atlas  and  Bartha  (4:393-394)  explain  that  a 
typical  petroleum  mixture  includes  aliphatics,  alicyclics,  aromatics  and  other  organics. 

The  microorganisms  that  metabolize  hydrocarbons  are  mostly  bacteria  and  fungi  (35 :308). 
Rainwater  and  Scholze  (43:108)  cite  a  study  by  Ghiorse  and  Balkwill  which  approximated 
the  number  of  bacteria  in  the  subsurface  at  one  million  per  gram  of  dry  soil.  Chapelle 
(18:336)  notes  that  Pseudomonas  putida  has  been  shown  in  studies  by  Gibson  and  others 
at  the  University  of  Texas  to  be  capable  of  degrading  benzene,  toluene  and  ethylbenzene. 
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Leahy  and  Colwell  (35:308)  listed  Archromobacter,  Acinetobacter,  Alcaligenes,  Arthrobacter, 
Bacillus,  Flavobacterium,  Nocardia  and  Pseudomonas  as  the  most  important  hydrocarbon 
degrading  bacterial  genera.  They  list  Trichoderma,  Mortierella,  Asperigillus  and  Penicillium 
as  the  most  important  fungal  genera  in  soil  media.  Other  genera  noted  in  the  literature  for 
their  hydrocarbon  degrading  ability  include  Xanthobacter  (54:1287)  and  Desulfobacterium 
(8:178).  Bak  and  Widel  isolated  a  new  species,  Desulfobacterium  phenolicolum,  an  obliga¬ 
tory  anaerobic  halophile,  capable  of  metabolizing  phenol.  Szewzyk  and  Pfennig  (53:164-165) 
identified  Desulfobacterium  catecholicum,  a  strictly  anaerobic,  sulfate  reducing  bacterium, 
capable  of  chemolithoautotrophic  growth  on  hydrogen  and  carbon  dioxide  or  heterotrophic 
growth  on  organics  such  as  catechol.  Heitkamp  and  Cemiglia  (29:1612)  isolated  a  halotoler- 
ant  strain  capable  of  mineralizing  naphthalene,  phenanthrene,  fluoranthrene,  pyrene  and  other 
aromatics  to  carbon  dioxide  when  supplied  with  nutrients  in  the  laboratory.  Cutright  and  Lee 
(22:404)  identified  a  species  of  Mycobacterium  capable  of  metabolizing  polycyclic  aromatic 
hydrocarbons. 

Although  protozoans  do  not  appear  to  utilize  hydrocarbons  directly  (35:309),  Madsen 
and  others  (37:252)  observed  high  density  protozoan  populations  growing  in  a  hydrocarbon 
plume.  The  protozoans  prey  on  bacteria  that  metabolize  the  hydrocarbons.  They  noted  that 
elevated  protozoan  biomass  acts  as  an  indicator  of  the  biodegradation  activity. 

After  a  petroleum  hydrocarbon  is  introduced  to  a  microbial  population  there  will  typi¬ 
cally  be  a  lag  time  before  biodegradation  begins  (58:1000)  (28:256).  During  this  period  the 
microbial  population  may  be  adapting  to  the  petroleum  hydrocarbon.  Leahy  and  Colwell 
(35:309)  define  adaptation  as  the  effect  of  prior  exposure  in  determining  the  rate  at  which 
a  microbial  community  will  biodegrade  hydrocarbons  and  cite  Spain  and  others  in  defining 
three  mechanisms  of  adaptation:  (1)  changes  in  the  action  of  enzymes,  (2)  mutations  to  genes 
allowing  new  metabolic  pathways  and  (3)  natural  selection  for  organisms  with  the  metabolic 
capability  to  degrade  the  petroleum  hydrocarbon.  They  note  that  DNA  encoded  on  plasmid 
may  play  an  important  role  in  genetic  adaptation  for  hydrocarbon  metabolism  because  of  its 
high  mobility  due  to  conjugation  and  transformation.  They  cited  a  study  by  Chakrabarty  on 
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Pseudomonas  that  showed  plasmid  DNA  was  encoded  for  metabolism  of  xylene,  toluene, 
naphthalene  and  other  compounds  suggesting  it  may  play  a  strong  role  in  adaptation. 

Zeyer  and  others  (65:946)  demonstrated  cross-acclimation,  defined  as  exposure  to  one 
compound  increasing  the  metabolism  of  other  compounds  of  similar  structure  (35:309),  by 
demonstrating  microorganisms  adapted  to  m-xylene  had  a  shorter  lag  time  than  expected  when 
metabolizing  toluene. 

Some  authors  report  no  adaptation.  Angley  and  others  (3:1406)  observed  no  lag  period 
in  lab  experiments  on  aerobic  degradation  of  BTEX  compounds  and  other  alkylbenzenes. 
Davis  and  others  (23:221)  observed  no  lag  time  in  experiments  on  aerobic  degradation  of 
benzene  until  initial  benzene  concentration  was  increased  from  about  1  ppm  to  10  ppm.  After 
the  increase  a  10  day  lag  occurred  before  biodegradation  began.  The  absence  of  a  lag  period 
indicates  the  microorganisms  have  no  need  to  adapt  to  the  substrate  before  they  can  begin 
metabolism. 

The  rate  at  which  metabolism  occurs  may  be  influenced  by  the  molecular  structure  of 
the  hydrocarbon.  Leahy  and  Colwell  (35:305)  cite  Perry  who  ranks  petroleum  compounds 
in  decreasing  order  of  susceptibility  to  biodegradation:  n-alkanes,  branched  alkanes,  low 
molecular  weight  aromatics,  and  cyclic  alkanes.  Chapelle  (18:350)  cites  a  study  by  Barker 
and  others  where  benzene,  toluene  and  xylenes  were  injected  into  an  aquifer.  They  found 
relative  degradation  rates  under  aerobic  conditions  to  be:  xylenes  >  toluene  >  benzene.  Under 
anaerobic  conditions  toluene  and  the  xylenes  were  almost  gone  by  108  days.  Benzene  was 
almost  gone  by  400  days. 

Stauffer  and  others  (5 1 :73)  supported  the  work  of  Barker  and  others  noted  above.  They 
injected  contaminants  into  groundwater  with  a  tritiated  water  tracer  and  monitored  contaminant 
concentrations.  They  verified  p-xylene  degrades  faster  than  benzene  under  aerobic  conditions 
(51:82) .  Their  results  showed  rapid  BTEX  biodegradation  under  aerobic  conditions. 

Angley  and  others  (3:1406)  performed  lab  experiments  on  aerobic  degradation  with 
BTEX  compounds  and  other  alkylbenzenes,  demonstrating  a  decrease  in  concentration  to 
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below  detection  limits  (a  decrease  of  four  orders  of  magnitude  for  the  BTEX  compounds) 
within  31  days.  This  was  accompanied  by  rapid  decreases  in  dissolved  oxygen  concentration. 

Davis  and  others  (23:221)  studied  biodegradation  in  the  lab.  They  showed  a  decrease 
in  benzene  concentration  of  50%  (from  about  1  ppm)  after  four  days  under  aerobic  conditions. 
The  biologically-inhibited  controls,  on  the  other  hand,  showed  only  2-18%  concentration 
decreases.  When  initial  benzene  was  increased  to  about  10  ppm,  50%  of  the  benzene  was 
gone  after  14  days  with  complete  degradation  by  35  days.  Under  sulfate  reducing  anaerobic 
conditions,  77  days  were  required  to  degrade  90%  of  1  ppm  benzene. 

Morgan  and  others  (40:189)  studied  both  aerobic  and  anaerobic  degradation  in  the 
lab.  They  found  that  under  aerobic  conditions  BTEX  rapidly  biodegraded.  Furthermore  they 
found  oxygen  supply  to  be  the  limiting  factor  for  degradation,  and  found  that  no  biodegradation 
occurred  under  anaerobic  conditions,  unless  nitrate  was  added  as  an  electron  acceptor. 

Baedecker  and  others  (6:569)  studied  samples  from  a  site  in  Minnesota  under  anaerobic 
conditions.  Their  data  showed  benzene  and  alkylbenzenes  degrade  under  aerobic  conditions. 
98%  degradation  of  benzene  and  toluene  occurred  in  125  and  45  days  respectively.  This 
was  coupled  with  increases  in  iron  and  manganese  concentrations,  suggesting  metal  reducing 
conditions. 

Neilsen  and  others  (41:461)  performed  a  laboratory  experiment  on  groundwater  and 
mixed  groundwater/sediment  media.  They  found  no  degradation  of  benzene,  toluene,  or 
o-xylene  under  anaerobic  conditions,  but  all  were  degraded  under  aerobic  conditions.  Degra¬ 
dation  of  benzene  and  toluene  required  40  to  70  days  and  resulted  in  reduction  to  less  than 
10%  of  initial  concentrations  (41:465). 

Other  authors  who  have  studied  this  subject  include  Wilson  and  others  (59:61)  whose 
lab  studies  confirmed  field  observations  that  aerobic  and  anaerobic  biodegradation  of  alkyl¬ 
benzenes  occurred.  Zeyer  and  others  (65:944)  showed  complete  mineralization  of  toluene  and 
m-xylene  under  anoxic  denitrifying  conditions  using  a  continuous  flow  system.  Wilson  and 
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others  (58:997)  demonstrated  methane  fermentation  (anaerobic  biodegradation)  of  toluene  in 
an  aquifer. 

The  process  of  biodegradation  is  one  of  enzymatic  transformation  through  a  series  of 
metabolic  intermediate  stages.  In  the  case  of  mineralization  the  end  products  are  carbon 
dioxide  and  water.  Kerfoot  (33:877)  gave  a  simplified  form  of  the  stoichiometric  equation  for 
aerobic  biodegradation: 


QH,  +  (a:  +  |)02  a;C02  +  (|)H20 

For  benzene,  the  degradation  process  involves,  first,  transformation  to  a  catechol  inter¬ 
mediate  and  then  breakdown  of  the  ring  structure  through  ortho  or  meta  cleavage.  Oxidation  of 
benzene  via  the  enzyme  benzene  dioxygenase  to  form  cis-benzene  dihydrodiol  is  the  first  step. 
The  second  step  involves  hydrogen  removal  via  the  enzyme  nicotine  adenine  dinucleotide 
(NAD).  A  diagram  of  this  process  is  shown  in  Figure  6. 


Catechol 

^OH 

^OH 

-  fn 

OH  /  ^ 

\  ^  OH 

V 

02 

NAD 

NADH2 

Benzene  dioxygenase) 

Figure  6.  Biodegradation  of  Benzene  to  Catechol 
(18:337) 
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The  ortho  cleavage  of  catechol  involves  breaking  the  ring  structure  between  the  two 
hydroxyl  groups.  Chapelle  (18:338)  describes  the  ortho  cleavage  pathway  through  several 
steps.  A  diagram  is  shown  in  Figure  7. 


Catechol 

nf 

Beta-ketoadipate-enol-lactone 

V 

^COOH 

(O'COOH 

SC^COOH 

^COOH 

^C=0 

(Further  degradation) 

Figure  7.  Ortho-Cleavage  Pathway 
(18:338) 

The  alternative  meta  cleavage  pathway  involves  ring  breakage  beside  one  of  the  two 
hydroxyl  groups  on  the  catechol,  followed  by  insertion  of  two  oxygen  atoms  via  the  enzyme 
catechol  2,3-dioxygenase.  The  ring  is  then  cleaved  by  a  hydrolase  (water  removing  enzyme). 
Further  degradation  of  the  organic  acid  product  to  CO2  may  then  occur.  A  diagram  is  shown 
in  Figure  8.  Note  that  beta-oxidation  in  the  figure  refers  to  oxidation  of  the  middle  carbon  and 
removal  as  carbon  dioxide. 

Chapelle  (18:337)  notes  that  ortho  and  meta  cleavage  of  catechol  has  been  identified 
in  Pseudomonas  putida,  Acinetobacter,  Bacillus,  Alcaligenes  and  Nocardia.  Breakdown  of 
toluene  and  xylenes  via  Pseudomonas  follow  processes  similar  to  that  of  benzene  except  that 
toluene  and  xylenes  form  methyl-  and  dimethylcatechol  respectively,  and  both  follow  the 
meta-cleavage  pathway  (18:339). 
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Figure  8.  Meta-Cleavage  Pathway 
(18:338) 


Anaerobic  transformation  of  benzene  and  toluene  was  studied,  using  gas  chromatog¬ 
raphy/mass  spectrometry  (28:254).  Results  suggest  initial  ring  hydroxylation  or  methyl- 
oxidation  (in  the  toluene  case)  leading  to  phenol,  cresol  and  aromatic  alcohol  intermediates. 
They  observed  incomplete  mineralization  to  methane  and  carbon  dioxide.  Cozzarelli  and  oth¬ 
ers  (20:138-140)  support  these  findings.  They  identified  oxygenated  intermediates  resulting 
from  anaerobic  degradation.  Benzene  produced  phenol,  toluene  produced  benzoic  acid  and 
xylenes  produced  toluic  acids.  The  theoretical  degradation  pathways  are  shown  in  Figures  9 
and  10. 

Cozzarelli  and  others  (21 :863)  found  that  benzene  and  alkylbenzenes  degrade  to  organic 
acid  intermediates  under  anaerobic  conditions.  They  also  found  that  reduction  of  iron,  man¬ 
ganese,  and  nitrogen  accompanied  oxidation  of  the  aromatics.  Anaerobic  biodegradation  by 
denitrification  is  insignificant  under  natural  conditions  due  to  the  usually  low  levels  of  nitrate 
in  most  contaminated  aquifers  (18:344). 

Benzene  degradation  under  methanogenic  conditions  begins  with  oxidation  of  benzene 
to  phenol,  through  several  aliphatic  acids  to  final  production  of  carbon  dioxide  and  methane 


55 


OH 

0 

\/ 

-"Aliphatic  acids 

Benzene  phenol 

Cyclohexane 

CO2 + CH^ 

Figure  9.  Anaerobic,  Ring  Hydroxylation  of  Benzene 
(28:258) 


(18:341).  Degradation  of  toluene  under  methanogenic  conditions  begins  with  oxidation  to 
benzylalcohol  and  then  to  benzaldehyde,  benzoic  acid,  and  through  several  aliphatic  acid  stages 
to  become  carbon  dioxide  and  methane.  Cozzarelli  and  others  (20:138-140)  may  support  this 
through  their  study  of  a  site  in  Minnesota.  Under  anaerobic  conditions,  they  monitored 
benzene  transformation  to  phenol  and  toluene  transformation  to  benzoic  acid.  Alternatively 
toluene  may  oxidize  directly  to  ortho  or  p-cresol  and  then  to  carbon  dioxide  and  methane 
(18:341). 

Under  iron  reducing  conditions  degradation  of  toluene  is  to  benzylalcohol,  benzalde¬ 
hyde,  benzoate,  and  finally  to  carbon  dioxide.  At  each  step  oxidation  of  the  aromatic  is 
accompanied  by  reduction  of  Fe(III)  to  Fe(II)  (18:343).  Baedecker  and  others  (6:580)  show 
a  theoretical  stoichiometric  equation  for  transformation  of  toluene  using  Iron  (III)  Hydroxide 
as  an  electron  acceptor: 

CtHs  +  36Fe(OH)3  -f  65H+  THCOg  -b  8Fe^  -f  87H2O 
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Figure  10.  Anaerobic,  Methyl-Oxidation  of  Toluene 
(28:258) 


The  oxygen  from  the  iron  (III)  hydroxide  was  inserted  into  the  ring  structure  of  toluene 
in  this  example.  In  other  anaerobic  cases  the  source  of  oxygen  inserted  into  ring  structures 
may  be  the  water  itself.  Vogel  and  Grbic-Galic  (55:200)  used  ^®0-labeled  water  to  show 
that  for  the  first  oxidation  step  in  the  anaerobic  transformation  of  benzene  and  toluene  under 
methanogenic  conditions,  the  oxygen  atom  came  from  water. 

Clearly,  a  variety  of  microbial  genera  are  capable  of  metabolizing  many  aromatic 
petroleum  hydrocarbon  compounds.  Many  conditions  once  thought  prohibitive  of  biodegra¬ 
dation  are  now  proven  possibilities  for  natural  attenuation.  Although  the  data  suggest  aerobic 
biodegradation  may  occur  more  rapidly  than  anaerobic  biodegradation,  it  is  apparent  that  both 
may  effectively  reduce  BTEX  concentrations  in  an  aquifer.  If  this  approach  received  more 
regulatory  acceptance,  savings  of  restoration  dollars  are  possible.  Remedial  costs  for  natural 
attenuation  candidate  sites  could  be  significantly  lower  than  under  alternative  technologies. 
If  natural  systems  can  be  self-cleaning,  public  funds  could  be  used  somewhere  other  than  on 
IRP  sites. 
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Appendix  B.  Initial  Bioplume  Input  File 

Biopluitie  Input  File  for  MADE-2  data 
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203.8 

203.9 

204.0 

204.1 

0.0 

0.0 

204.2 

204.2 

204.2 

204.2 

204.2 

204.2 

204.2 

204.2 

204.2 

204.1 

204.1 

204.0 

203.9 

203.9 

203.9 

204.0 

204.1 

204.1 

0.0 

0.0 

204.2 

204.2 

204.2 

204.2 

204.3 

204.3 

204.3 

204.3 

204.3 

204.3 

204.2 

204.2 

204.1 

204.0 

204.1 

204.1 

204.2 

204.2 

0.0 

0.0 

204.3 

204.3 

204.3 

204.3 

204.3 

204.4 

204.4 

204.4 

204.4 

204.4 

204.4 

204.3 

204.2 

204.2 

204.2 

204.2 

204.3 

204.3 

0.0 

0.0 

204.3 

204.3 

204.4 

204.4 

204.4 

204.4 

204.5 

204.5 

204.5 

204.6 

204.5 

204.4 

204.4 

204.4 

204.4 

204.4 

204.4 

204.4 

0.0 

0.0 

204.4 

204.4 

204.4 

204.5 

204.5 

204.5 

204.6 

204.6 

204.8 

204.8 

204.6 

204.6 

204.5 

204.5 

204.5 

204.5 

204.5 

204.5 

0.0 

0.0 

204.4 

204.5 

204.5 

204.5 

204.6 

204.6 

204.7 

204.7 

204.8 

204.8 

204.7 

204.7 

204.6 

204.6 

204.6 

204.6 

204.6 

204.6 

0.0 

0.0 

204.5 

204.5 

204.5 

204.6 

204.6 

204.7 

204.7 

204.8 

204.8 

204.9 

204.8 

204.8 

204.7 

204.7 

204.8 

204.8 

204.7 

204.7 

0.0 

0.0 

0  0.0 

0.0 

0.0 

0.0 

0.0 

0.0 

0.0 

0.0 

0.0 

0.0 

0.0  ' 

0.0  ( 

).0  C 

1.0  0 

.0  0 

.0  0, 

.0  0. 

0  0. 

0  0.0 

1  1000. 


0.00 

0.00 

0.00 

0.00 

0.00 

0.00 

0.00 

0.00 

0.00 

0.00 

0.00 

0.00 

0.00 

0.00 

0.00 

0.00 

0.00 

0.00 

0.00 

0.00 

0.00 

5.37 

5.41 

5.44 

5.48 

5.50 

5.52 

5.53 

5.50 

5.42 

5.24 

4.82 

4.12 

3.36 

3.11 

3.16 

3.24 

3.31 

3.38 

0.00 

0.00 

5.35 

5.39 

5.43 

5.47 

5.50 

5.52 

5.53 

5.50 

5.43 

5.27 

4.85 

4.08 

3.30 

3.06 

3.14 

3.22 

3.31 

3.38 

0.00 

0.00 

5.33 

5.37 

5.41 

5.45 

5.48 

5.51 

5.52 

5.49 

5.41 

5.26 

4.92 

4.07 

3.44 

3.18 

3.19 

3.26 

3.34 

3.41 

0.00 

0.00 

5.30 

5.34 

5.38 

5.42 

5.45 

5.48 

5.49 

5.47 

5.37 

5.12 

4.62 

3.98 

3.54 

3.32 

3.29 

3.33 

3.40 

3.46 

0.00 

0.00 

5.26 

5.30 

5.34 

5.38 

5.41 

5.44 

5.46 

5.44 

5.33 

4.96 

4.20 

3.84 

3.61 

3.46 

3.41 

3.44 

3.48 

3.53 

0.00 

0.00 

5.22 

5.26 

5.29 

5.33 

5.36 

5.39 

5.41 

5.40 

5.34 

5.00 

4.13 

3.85 

3.70 

3.59 

3.55 

3.56 

3.59 

3.62 

0.00 

0.00 

5.17 

5.20 

5.24 

5.27 

5.30 

5.32 

5.33 

5.33 

5.32 

5.17 

4.39 

4.00 

3.83 

3.74 

3.70 

3.69 

3.71 

3.73 

0.00 

0.00 

5.11 

5.15 

5.18 

5.20 

5.23 

5.24 

5.24 

5.22 

5.17 

4.96 

4.49 

4.15 

3.98 

3.89 

3.85 

3.84 

3.83 

3.84 

0.00 

0.00 

5.06 

5.08 

5.11 

5.13 

5.14 

5.15 

5.13 

5.08 

4.98 

4.78 

4.49 

4.26 

4.12 

4.05 

4.01 

3.99 

3.97 

3.96 

0.00 

0.00 

4.99 

5.01 

5.04 

5.05 

5.06 

5.05 

5.02 

4.95 

4.81 

4.61  4.46 

4.33 

4.25 

4.20 

4.17 

4.14 

4.11 

4.08 

0.00 

0.00 

4.92 

4.94 

4.96 

4.97 

4.97 

4.96 

4.92 

4.83 

4.68 

4.46 

4.44 

4.41 

4.38 

4.36 

4.33 

4.30 

4.25 

4.20 

0.00 

0.00 

4.85 

4.87 

4.88 

4.89 

4.89 

4.87 

4.83 

4.75 

4.62 

4.49 

4.48 

4.50 

4.51 

4.51 

4.49 

4.46 

4.40 

4.33 

0.00 

0.00 

4.78 

4.79 

4.80 

4.80 

4.80 

4.79 

4.76 

4.70 

4.61 

4.54 

4.55 

4.60 

4.64 

4.67 

4.66 

4.61 

4.53 

4.44 

0.00 

0.00 

4.70 

4.71 

4.72 

4.72 

4.72 

4.72 

4.70 

4.67 

4.62 

4.59 

4.64 

4.71 

4.78 

4.82 

4.83 

4.77 

4.66 

4.54 

0.00 

0.00 

4.62 

4.63 

4.63 

4.63 

4.64 

4.64 

4.64 

4.64 

4.64 

4.62 

4.73 

4.82 

4.89 

4.96 

4.99 

4.91 

4.77 

4.63 

0.00 

0.00 

4.54 

4.54 

4.54 

4.55 

4.55 

4.57 

4.59 

4.62 

4.67 

4.74 

4.86 

4.91 

4.98 

5.08 

5.15 

5.03 

4.85 

4.68 

0.00 

0.00 

4.46 

4.45 

4.45 

4.45 

4.46 

4.48 

4.52 

4.59 

4.69 

4.85 

4.97 

4.95 

5.00 

5.12 

5.28 

5.08 

4.88 

4.70 

0.00 

0.00 

4.38 

4.36 

4.36 

4.35 

4.36 

4.38 

4.43 

4.52 

4.66 

4.83 

4.92 

4.90 

4.95 

5.05 

5.13 

5.02 

4.85 

4.69 

0.00 

0.00 

4.29 

4.27 

4.26 

4.25 

4.25 

4.27 

4.32 

4.42 

4.57 

4.70 

4.73 

4.76 

4.82 

4.90 

4.95 

4.90 

4.78 

4.64 

0.00 

0.00 

4.20 

4.18 

4.16 

4.14 

4.13 

4.13 

4.17 

4.26 

4.42 

4.56 

4.54 

4.56 

4.63 

4.72 

4.77 

4.75 

4.67 

4.57 

0.00 

0.00 

4.12 

4.09 

4.06 

4.03 

4.00 

3.98 

3.99 

4.03 

4.19 

4.41 

4.25 

4.31 

4.41 

4.51 

4.58 

4.59 

4.55 

4.48 

0.00 

0.00 

4.04 

4.00 

3.96 

3.91 

3.87 

3.82 

3.79 

3.76 

3.77 

3.71 

3.85 

4.02 

4.18 

4.31 

4.39 

4.43 

4.42 

4.38 

0.00 

0.00 

3.96 

3.91 

3.86 

3.80 

3.74 

3.67 

3.59 

3.48 

3.31 

3.19 

3.49 

3.76 

3.96 

4.12 

4.22 

4.28 

4.29 

4.28 

0.00 

0.00 

3.88 

3.83 

3.77 

3.70 

3.63 

3.54 

3.42 

3.27 

3.05 

2.96 

3.27 

3.55 

3.77 

3.94 

4.06 

4.13 

4.17 

4.18 

0.00 

0.00 

3.81 

3.75 

3.69 

3.61 

3.53 

3.43 

3.30 

3.15 

2.99 

2.97 

3.17 

3.41 

3.62 

3.79 

3.92 

4.00 

4.05 

4.08 

0.00 

0.00 

3.75 

3.68 

3.61 

3.54 

3.45 

3.34 

3.22 

3.09 

2.99 

2.99 

3.13 

3.32 

3.51 

3.67 

3.80 

3.89 

3.95 

3.98 

0.00 

0.00 

3.69 

3.62 

3.55 

3.47 

3.38 

3.28 

3.18 

3.07 

3.00 

3.01 

3.11 

3.27 

3.43 

3.58 

3.70 

3.79 

3.85 

3.90 

0.00 

0.00 

3.63 

3.57 

3.49 

3.42 

3.33 

3.24 

3.15 

3.07 

3.02 

3.03 

3.11 

3.23 

3.37 

3.50 

3.61 

3.70 

3.77 

3.82 

0.00 

0.00 

0.00 

0.00 

0.00 

0.00 

0.00 

0.00 

0.00 

0.00 

0.00 

0.00 

0.00 

0.00 

0.00 

0.00 

0.00 

0.00 

0.00 

0.00 

0.00 

1 

92  92  7  200  1  0  0  0  0  0  0.2521  0  0 

11  23  0.0  0.0  0.0 

0 

0 

0 

0 
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Appendix  C.  SAS  Output 

C.  1  Stepwise  Procedure  for  RSM  Screening  Phase 

SAS  Program: 

option  linesize=80 ; 
data  dl; 

infile  ' rsmphsl . dat ' ; 
input  t  p  d  y; 
tp=t*p; 
td=t*d; 
pd=p*d; 
tpd=t*p*d; 
proc  print; 
proc  stepwise; 

model  y  =  t  p  d  tp  td  pd  tpd/stepwise; 
run; 


SAS  Output: 


The  SAS  System 

08:51  Friday,  October  6 


OBS 

T 

P 

D 

Y 

TP 

TD 

PD 

TPD 

1 

-1 

-1 

-1 

27246.50 

1 

1 

1 

-1 

2 

-1 

-1 

1 

544968.00 

1 

-1 

-1 

1 

3 

-1 

1 

-1 

27246.50 

-1 

1 

-1 

1 

4 

-1 

1 

1 

544968.00 

-1 

-1 

1 

-1 

5 

1 

-1 

-1 

2.22 

-1 

-1 

1 

1 

6 

1 

-1 

1 

1.81 

-1 

1 

-1 

-1 

7 

1 

1 

-1 

2.22 

1 

-1 

-1 

-1 

8 

1 

1 

1 

1.81 

1 

1 

1 

1 

The  SAS  System 

08:51  Friday,  October  6 
Stepwise  Procedure  for  Dependent  Variable  Y 


16 

,  1995 


17 

,  1995 
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step  1  Variable  T  Entered  R-square  =  0.37918514  C(p)  = 

DF  Sum  of  Squares  Mean  Square  F  Prob>F 

Regression  1  163712412173.29  163712412173.29  3.66  0.1041 

Error  6  268035551562.41  44672591927.069 

Total  7  431747963735.70 

Parameter  Standard  Type  II 

Variable  Estimate  Error  Sum  of  Squares  F  Prob>F 

INTERCEP  143054.63198250  74726.66184759  163717021853.19  3.66  0.1041 

T  -143052.6180175  74726.66184759  163712412173.29  3.66  0.1041 

Bounds  on  condition  number:  1,  1 


The  SAS  System  18 

08:51  Friday,  October  6,  1995 

Step  2  Variable  TD  Entered  R-square  =  0.68959305  C(p)  = 

DF  Sum  of  Squares  Mean  Square  F  Prob>F 

Regression  2  297730396829.23  148865198414.62  5.55  0.0537 

Error  5  134017566906.47  26803513381.293 

Total  7  431747963735.70 

Parameter  Standard  Type  II 

Variable  Estimate  Error  Sum  of  Squares  F  Prob>F 

INTERCEP  143054.63198250  57882.97826358  163717021853.19  6.11  0.0564 

T  -143052.6180175  57882.97826358  163712412173.29  6.11  0.0564 

TD  -129430.4758625  57882.97826358  134017984655.95  5.00  0.0756 

Bounds  on  condition  number:  1,  4 

The  SAS  System  19 

08:51  Friday,  October  6,  1995 

Step  3  Variable  D  Entered  R-square  =  1.00000000  C(p)  = 

DF  Sum  of  Squares  Mean  Square  F  Prob.>F 
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Regression  3  431747963735.70  143915987911.90 

Error  4  0.00000000  0.00000000 

Total  7  431747963735.70 

Parameter  Standard  Type  II 

Variable  Estimate  Error  Sum  of  Squares  F  Prob>F 

INTERCEP  143054.63198250  0.00000000  163717021853.19 

T  -143052.6180175  0.00000000  163712412173.29 

D  129430.27413750  0.00000000  134017566906.47 

TD  -129430.4758625  0.00000000  134017984655.95 

Bounds  on  condition  number:  1,  9 


All  variables  left  in  the  model  are  significant  at  the  0.1500  level. 

No  other  variable  met  the  0.1500  significance  level  for  entry  into  the  model. 

The  SAS  System  20 

08:51  Friday,  October  6,  1995 

Summary  of  Stepwise  Procedure  for  Dependent  Variable  Y 

Variable  Number  Partial  Model 

Step  Entered  Removed  In  R**2  R**2  C(p)  F  Prob>F 

IT  1  0.3792  0.3792  .  3.6647  0.1041 

2  TD  2  0.3104  0.6896  .  5.0000  0.0756 

3D  3  0.3104  1.0000 

C.  2  First  Run  of  the  First-  Order  Design  Phase 

SAS  Program: 

option  linesize=80 ; 
data  dl ; 

infile  '  rsmphs2  .  dat '  ; 
input  t  d  y; 
proc  print; 
proc  reg; 

model  y  =  t  d; 
run; 
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SAS  Output: 


The  SAS  System  7 

08:51  Friday,  October  6,  1995 


OBS 

T 

D 

Y 

1 

1 

1 

2.11605 

2 

1 

-1 

2.36023 

3 

-1 

1 

1.24133 

4 

-1 

-1 

3.37311 

The 

SAS  System 

08:51  Friday,  October  6,  1995 

Model:  MODEL 1 
Dependent  Variable:  Y 

Analysis  of  Variance 


Source 

Sum  of 

DF  Squares 

Mean 

Square 

F  Value 

Prob>F 

Model 

Error 

C  Total 

2  1.41607 

1  0.89076 

3  2.30683 

0.70803 

0.89076 

0.795 

0.6214 

Root  MSE 

Dep  Mean 

C.V. 

0.94380  R-square 

2.27268  Adj  R-sq 

41.52806 

0.6139 

-0.1584 

The  SAS  System  9 

08:51  Friday,  October  6,  1995 

Parameter  Estimates 

Parameter  Standard  T  for  HO : 

Variable  DF  Estimate  Error  Parameter=0  Prob  >  |t| 
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INTERCEP  1  2.272680  0.47190000  4.816  0.1303 
T  1  -0.034540  0.47190000  -0.073  0.9535 
D  1  -0.593990  0.47190000  -1.259  0.4274 


C.3  40%  Design  without  Centerpoint 
SAS  Program: 

option  linesize=80; 
data  dl; 

infile  '  rsmphs3  .  dat '  ; 
input  t  d  y; 
proc  print; 
proc  reg; 

model  y  =  t  d; 
run; 

SAS  Output: 

The  SAS  System  1 

17:11  Thursday,  October  26,  1995 


OBS 

T 

D 

Y 

1 

-1 

-1 

2.21077 

2 

-1 

1 

7.46751 

3 

1 

-1 

0.96785 

4 

1 

1 

2.21034 

The  SAS  System  2 

17:11  Thursday,  October  26,  1995 

Model:  MODEL 1 
Dependent  Variable:  Y 

Analysis  of  Variance 
Sum  of  Mean 

Source  DF  Squares  Square  F  Value  Prob>F 

Model  2  21.12278  10.56139  2.622  0.4002 

Error  1  4.02855  4.02855 
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C  Total 


3 


25.15133 


Root  MSE 
Dep  Mean 

C.V. 


2.00713 

3.21412 

62.44718 


R-square 
Adj  R-sq 


0.8398 

0.5195 


The  SAS  System 


3 


17:11  Thursday,  October  26,  1995 


Parameter  Estimates 


Variable  DF 

INTERCEP  1 
T  1 
D  1 


Parameter 

Estimate 

3.214118 

-1.625022 

1.624807 


Standard 

Error 

1.00356300 

1.00356300 

1.00356300 


T  for  HO: 
Parameter=0 

3.203 

-1.619 

1.619 


Prob  >  I T I 

0.1927 

0.3522 

0.3522 


C.4  40%  Design  With  Centerpoint 

SAS  Program: 

option  linesize=80; 
data  dl; 

infile  'rsmphs3.dat'; 
input  t  d  y; 
proc  print; 
proc  reg; 

model  y  =  t  d; 
run; 

SAS  Output: 


The  SAS  System  21 

08:51  Friday,  October  6,  1995 

OBS  T  D  Y 

1  -1  -1  2.21077 

2-11  7.46751 
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3 

4 

5 


1  -1  0.96785 
1  1  2.21034 
0  0  2.21077 

The  SAS  System  22 
08:51  Friday,  October  6,  1995 


Model:  MODELl 
Dependent  Variable:  Y 


Analysis  of  Variance 


Sum  of 

Mean 

Source 

DF 

Squares 

Square 

F  Value 

Prob>F 

Model 

2 

21.12278 

10.56139 

4.370 

0.1862 

Error 

2 

4.83392 

2.41696 

C  Total 

4 

25.95670 

Root  MSB  1.55466  R-square  0.8138 
Dep  Mean  3.01345  Adj  R-sq  0.6275 
C.V.  51.59065 


The  SAS  System 


23 

08:51  Friday,  October  6,  1995 


Parameter  Estimates 


Variable 

DF 

Parameter 

Estimate 

Standard 

Error 

T  for  HO: 

Parameter=0 

Prob  >  |t| 

INTERCEP 

1 

3.013448 

0.69526402 

4.334 

0.0493 

T 

1 

-1.625022 

0.77732880 

-2.091 

0.1717 

D 

1  ' 

1.624807 

0.77732880 

2.090 

0.1718 

C.5  5%  Design  Around  50%  Step  along  Gradient 
SAS  Program : 

option  linesize=80; 
data  dl; 

infile  ' Bdesign . dat ' ; 
input  t  d  y; 
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proc  print; 
proc  reg; 

model  y  =  t  d; 
run; 

SAS  Output : 

The  SAS  System  16 

15:55  Saturday,  October  7,  1995 


OBS 

T 

D 

Y 

1 

-1 

-1 

0.94093 

2 

-1 

1 

1.01783 

3 

1 

-1 

0.91888 

4 

1 

1 

0.94280 

5 

0 

0 

0.94190 

The  SAS  System  17 

15:55  Saturday,  October  7,  1995 

Model:  MODEL 1 
Dependent  Variable:  Y 

Analysis  of  Variance 


Sum  of 

Mean 

Source 

DF  Squares 

Square 

F  Value 

Prob>F 

Model 

Error 

C  Total 

2  0.00490 

2  0.00084 

4  0.00574 

0.00245 

0.00042 

5.821 

0.1466 

Root  MSE 

0.02051  R-square 

0.8534 

Dep  Mean 

C.V. 

0.95247  Adj 

2.15340 

R-sq 

0.7068 

The  SAS  System  18 

15:55  Saturday,  October  7,  1995 


Parameter  Estimates 
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Parameter  Standard  T  for  HO : 

Variable  DF  Estimate  Error  Parameter=0  Prob  >  |t| 

INTERCEP  1  0.952465  0.00917254  103.839  0.0001 

T  1  -0.024272  0.01025521  -2.367  0.1416 

D  1  0.025206  0.01025521  2.458  0.1332 

C.6  5%  Design  on  Step  9  Point 
SAS  Program: 

option  linesize=80; 
data  dl ; 

infile  '  Bdesign .  dat '  ; 
input  t  d  y; 
proc  print; 
proc  reg; 

model  y  =  t  d; 
run; 

SAS  Output: 


The  SAS  System  1 

15:37  Sunday,  October  8,  1995 


OBS 

T 

D 

Y 

1 

-1 

-1 

0.92522 

2 

-1 

1 

0.97780 

3 

1 

-1 

0.92392 

4 

1 

1 

0.92522 

5 

0 

0 

0.92522 

The  SAS  System  2 

15:37  Sunday,  October  8,  1995 

Model :  MODEL 1 
Dependent  Variable:  Y 

Analysis  of  Variance 
Sum  of  Mean 

Source  DF  Squares  Square  F  Value  Prob>F 


Model  2  0.00145  0.00073  1.840  0.3521 

Error  2  0.00079  0.00039 

C  Total  4  0.00224 

Root  MSE  0.01986  R-square  0.6479 

Dep  Mean  0.93547  Adj  R-sq  0.2959 

C.V.  2.12296 

The  SAS  System  3 

15:37  Sunday,  October  8,  1995 

Parameter  Estimates 

Parameter  Standard  T  for  HO : 

Variable  DF  Estimate  Error  Parameter=0  Prob  >  |t| 

INTERCEP  1  0.935474  0.00888155  105.328  0.0001 

T  1  -0.013471  0.00992988  -1.357  0.3077 

D  1  0.013471  0.00992988  1.357  0.3077 

C.7  First  Ridge  Analysis 
SAS  Program: 

linesize=80 ; 
data  dl; 

infile='phs2 .dat' ; 
input  t  d  y; 
proc  print; 
proc  rsreg; 

model  y  =  t  d/nocode; 
ridge  min  radius  =  0  to  4  by  0.1; 
run; 


SAS  Output: 


The  SAS  System  16:34  Sunday,  October  8,  1995  4 
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OBS 

T 

D 

Y 

TD 

TT 

DD 

1 

-1.000 

-1.000 

0.92522 

1 

1.00000 

1.00000 

2 

-1.000 

1.000 

0.97780 

-1 

1.00000 

1.00000 

3 

1.000 

-1.000 

0.92392 

-1 

1.00000 

1.00000 

4 

1.000 

1.000 

0.92522 

1 

1.00000 

1.00000 

5 

0.000 

0.000 

0.92522 

0 

0.00000 

0.00000 

6 

1.414 

0.000 

0.91949 

0 

1.99940 

0.00000 

7 

-1.414 

0.000 

0.95819 

0 

1.99940 

0.00000 

8 

0.000 

1.414 

0.95497 

0 

0.00000 

1.99940 

9 

0.000 

-1.414 

0.91991 

0 

0.00000 

1.99940 

The  S AS  System  16:34  Sunday,  October  8,  1995  E 


Response  Surface  for  Variable  Y 


Response  Mean 
Root  MSE 
R-Square 

Coef.  of  Variation 


0.936658 

0.000898 

0.9993 

0.0959 


Regression 


Degrees 

of 

Type  I  Sum 

Freedom 

of  Squares 

R- Square 

F-Ratio 

Prob  >  F 

Linear 
Quadratic 
Crossproduct 
Total  Regress 


2 

2 

1 

5 


0.002813 
0.000149 
0.000657 
0.003619 
The  SAS  System 


0.7766 

0.0412 

0.1815 

0.9993 


1744.2 

92.520 

815.3 

897.8 


0.0000 

0.0020 

0.0001 

0.0001 


16:34  Sunday,  October  8,  1995 


e 


Residual 


Degrees 

of  Sum  of 

Freedom  Squares  Mean  Square 


Total  Error 


0.000002419  0.000000806 


Degrees 

of  Parameter 


Standard  T  for  HO : 
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Parameter 

Freedom 

Estimate 

Error 

Parameter 

-0  Prob  > 

INTERCEPT 

1 

0.925218 

0.000898 

1030.4 

0.0000 

T 

1 

-0.013577 

0.000317 

-42.765 

0.0000 

D 

1 

0.012934 

0.000317 

40.739 

0.0000 

1 

0.006785 

0.000527 

12.887 

0.0010 

D*T 

1 

-0.012819 

0.000449 

-28.553 

0.0001 

D*D 

1 

0.006087 

0.000527 

11.560 

0.0014 

The  SAS  System  16:34 

Sunday,  October  8,  1995 

Degrees 

of 

Sum  of 

Factor 

Freedom 

Squares 

Mean  Square 

F-Ratio 

Prob  >  F 

T 

3 

0.002266 

0.000755 

936.7 

0.0001 

D 

3 

0.002103 

0.000701 

869.5 

0.0001 

The  SAS  System  16:34  Sunday,  October  8,  1995  E 


Canonical  Analysis  of  Response  Surface 


Critical 
Factor  Value 

T  -0.595868 

D  -1.689911 


Predicted  value  at  stationary  point  0.918335 


Eigenvalues 


Eigenvectors 
T  D 


0.012855 

0.000016930 


0.726087  -0.687603 

0.687603  0.726087 


Stationary  point  is  a  minimum. 

The  SAS  System  16:34  Sunday,  October  8,  1995  S 
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Estimated  Ridge  of  Minimum  Response  for  Variable  Y 


Estimated 

Standard 

Factor  Values 

Radius 

Response 

Error 

T 

D 

0 

0.925218 

0.000898 

0 

0 

0.100000 

0.923471 

0.000894 

0.072373 

-0.069009 

0.200000 

0.921982 

0.000882 

0.144656 

-0.138110 

0.300000 

0.920750 

0.000863 

0.216789 

-0.207371 

0.400000 

0.919774 

0.000837 

0.288631 

-0.276933 

0.500000 

0.919056 

0.000804 

0.359806 

-0.347188 

0.600000 

0.918595 

0.000767 

0.428764 

-0.419716 

0.700000 

0.918390 

0.000726 

0.474604 

-0.514540 

0.800000 

0.918363 

0.000684 

0.295297 

-0.743505 

0.900000 

0.918355 

0.000645 

0.163626 

-0.885001 

1.000000 

0.918350 

0.000615 

0.057195 

-0.998363 

1.100000 

0.918346 

0.000599 

-0.037907 

-1.099347 

1.200000 

0.918343 

0.000604 

-0.126535 

-1.193310 

1.300000 

0.918340 

0.000637 

-0.210988 

-1.282764 

1.400000 

0.918338 

0.000699 

-0.292543 

-1.369094 

1.500000 

0.918336 

0.000789 

-0.371980 

-1.453145 

1.600000 

0.918335 

0.000904 

-0.449810 

-1.535471 

The 

SAS  System 

16:34  Sunday, 

October  8,  1995 

Estimated 

Standard 

Factor  Values 

Radius 

Response 

Error 

T 

D 

1.700000 

0.918335 

0.001040 

-0.526389 

-1.616451 

1.800000 

0.918335 

0.001196 

-0.601967 

-1.696359 

1.900000 

0.918335 

0.001368 

-0.676733 

-1.775396 

2.000000 

0.918335 

0.001556 

-0.750829 

-1.853714 

2.100000 

0.918337 

0.001757 

-0.824363 

-1.931431 

2.200000 

0.918338 

0.001971 

-0.897423 

-2.008639 

2.300000 

0.918340 

0.002198 

-0.970077 

-2.085414 

2.400000 

0.918342 

0.002437 

-1.042383 

-2.161814 

2.500000 

0.918344 

0.002688 

-1.114384 

-2.237889 

2.600000 

0.918347 

0.002950 

-1.186121 

-2.313681 

2.700000 

0.918350 

0.003224 

-1.257624 

-2.389222 

2.800000 

0.918354 

0.003509 

-1.328920 

-2.464543 

2.900000 

0.918358 

0.003805 

-1.400033 

-2.539667 

3.000000 

0.918362 

0.004112 

-1.470982 

-2.614615 
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3.100000 

3.200000 

3.300000 

3.400000 

3.500000 

3.600000 


0.918367  0.004430  -1.541783  -2.689406 
0.918372  0.004759  -1.612451  -2.764055 
0.918377  0.005099  -1.683000  -2.838576 
0.918383  0.005449  -1.753440  -2.912979 
0.918389  0.005811  -1.823780  -2.987277 
0.918395  0.006183  -1.894030  -3.061478 

The  SAS  System  16:34  Sunday,  October  8,  1995  1 


Estimated  Standard 

Radius  Response  Error 


Factor  Values 
T  D 


3.700000 

0.918402 

0 

3.800000 

0.918409 

0 

3.900000 

0.918416 

0 

4.000000 

0.918424 

0 

006565  -1.964197  -3.135591 
006959  -2.034288  -3.209622 
007363  -2.104309  -3.283578 
007777  -2.174265  -3.357465 


C.8  Screening  Design  for  Transport  Calibration 

SAS  Program: 

option  linesi2e=80; 
data  dl; 

infile  'transl.dat'; 
input  alphaL  ratio  F  Rf  RC  y; 
proc  print; 
proc  reg; 

model  y  =  alphaL  ratio  F  Rf  RC; 
run; 

SAS  Output : 

The  SAS  System  1 

08:52  Wednesday,  October  18,  1995 


OBS 

ALPHAL 

RATIO 

F 

RF 

RC 

Y 

1 

-1 

-1 

-1 

-1 

-1 

94.047 

2 

-1 

-1 

-1 

-1 

1 

67.183 

3 

-1 

-1 

-1 

1 

-1 

219.256 

4 

-1 

-1 

-1 

1 

1 

77.893 

5 

-1 

-1 

1 

-1 

-1 

172.274 

6 

-1 

-1 

1 

-1 

1 

73.619 
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7 

-1 

-1 

1 

1 

-1 

273.971 

8 

-1 

-1 

1 

1 

1 

93.629 

9 

-1 

1 

-1 

-1 

-1 

67.183 

10 

-1 

1 

-1 

-1 

1 

67.183 

11 

-1 

1 

-1 

1 

-1 

67.183 

12 

-1 

1 

-1 

1 

1 

67.183 

13 

-1 

1 

1 

-1 

-1 

67.183 

14 

-1 

1 

1 

-1 

1 

67.183 

15 

-1 

1 

1 

1 

-1 

67.174 

16 

-1 

1 

1 

1 

1 

67.183 

17 

1 

-1 

-1 

-1 

.  -1 

67.183 

18 

1 

-1 

-1 

-1 

1 

67.183 

19 

1 

-1 

-1 

1 

-1 

78.544 

20 

1 

-1 

-1 

1 

1 

67.183 

The 

SAS  System 

2 

08: 

52  Wednesday,  October  18, 

1995 

OBS 

ALPHAL 

RATIO  F 

RF 

RC 

Y 

21 

1 

-1 

1 

-1 

-1 

67.051 

22 

1 

-1 

1 

-1 

1 

67.183 

23 

1 

-1 

1 

1 

-1 

128.365 

24 

1 

-1 

1 

1 

1 

67.949 

25 

1 

1 

-1 

-1 

-1 

67.183 

26 

1 

1 

-1 

-1 

1 

67.183 

27 

1 

1 

-1 

1 

-1 

67.183 

28 

1 

1 

-1 

1 

1 

67.183 

29 

1 

1 

1 

-1 

-1 

67.183 

30 

1 

1 

1 

-1 

1 

67.183 

31 

1 

1 

1 

1 

-1 

67.183 

32 

1 

1 

1 

1 

1 

67.183 

The 

SAS  System 

3 

08: 

52  Wednesday,  October  18, 

1995 

Model:  MODELl 

Dependent  Variable 

:  Y 

Analysis  of  Variance 
Sum  of  Mean 

Source  DF  Squares  Square  F  Value  Prob>F 
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Model 
Error 
C  Total 


5  31342.60698  6268.52140 
26  39503.12918  1519.35112 
31  70845.73616 


4.126 


0.0068 


Root  MSE  38.97885  R-square  0.4424 

Dep  Mean  86.16980  Adj  R-sq  0.3352 

C.V.  45.23494 

The  SAS  System  4 

08:52  Wednesday,  October  18,  1995 

Parameter  Estimates 

Parameter  Standard  T  for  HO : 


Variable 

DF 

Estimate 

Error 

Parameter=0 

Prob  >  |t| 

INTERCEP 

1 

86.169800 

6.89055314 

12.505 

0.0001 

ALPHAL 

1 

-14.413200 

6.89055314 

-2.092 

0.0464 

RATIO 

1 

-18.987263 

6.89055314 

-2.756 

0.0106 

F 

1 

6.423706 

6.89055314 

0.932 

0.3598 

RF 

1 

10.345581 

6.89055314 

1.501 

0.1453 

RC 

1 

-16.214344 

6.89055314 

-2.353 

0.0265 

C.9  First-Order  Design  for  RSM  Transport  Calibration 
SAS  Program: 

option  linesize=80; 
data  dl; 

infile  'trans2.dat'; 
input  alphaL  ratio  y; 
y=y-67; 
proc  print; 
proc  reg; 

model  y  =  alphaL  ratio; 
run; 


SAS  Output: 

The  SAS  System  4 

17:24  Wednesday,  October  18,  1995 
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OBS 

ALPHAL 

RATIO 

Y 

1 

-1 

-1 

0.1906 

2 

-1 

1 

0.0921 

3 

1 

-1 

0.0505 

4 

1 

1 

0.1115 

The  SAS  System  5 

17:24  Wednesday,  October  18,  1995 


Model:  MODELl 
Dependent  Variable:  Y 


Source  DF 

Model  2 
Error  1 
C  Total  3 


Analysis  of  Variance 


Sum  of 

Mean 

Squares 

Square 

0.00399 

0.00200 

0.00636 

0.00636 

0.01035 

F  Value  Prob>F 

0.314  0.7838 


Root  MSE 
Dep  Mean 

C.V. 


0.07975 

0.11118 

71.73375 


R-square 
Adj  R-sq 


0.3857 

-0.8428 


The  SAS  System  6 

17:24  Wednesday,  October  18,  1995 

Parameter  Estimates 


Variable 

DF 

Parameter 

Estimate 

Standard 

Error 

T  for  HO : 

Parameter=0 

Prob  >  |t| 

INTERCEP 

1 

0.111175 

0.03987500 

2.788 

0.2192 

ALPHAL 

1 

-0.030175 

0.03987500 

-0.757 

0.5876 

RATIO 

1 

-0.009375 

0.03987500 

-0.235 

0.8530 

C.IO  Ridge  Analysis  for  RSM  Transport  Calibration 
SAS  Program: 

option  linesize=80 ; 
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data  dl; 

infile  'trans3.dat'; 
input  alphaL  ratio  y; 
proc  print; 
proc  rsreg; 

model  y  =  alphaL  ratio/nocode; 
ridge  min  radius  =  0  to  5  by  0 . 1  ; 
run; 


SAS  Output : 

The  SAS  System  7 

17:11  Thursday,  October  26,  1995 


OBS 

ALPHAL 

RATIO 

Y 

1 

-1.000 

-1.000 

67.1906 

2 

-1.000 

1.000 

67.0921 

3 

1.000 

-1.000 

67.0505 

4 

1.000 

1.000 

67.1115 

5 

-1.414 

0.000 

67.2464 

6 

1.000 

0.000 

67.0764 

7 

0.000 

1.414 

67.0327 

8 

0.000 

-1.000 

67.0394 

9 

0.000 

0.000 

67.0362 

The 

SAS  System 

17:11  Thursday,  October  26,  1995 


Response  Surface  for  Variable  Y 


Response  Mean 
Root  MSE 
R- Square 

Coef.  of  Variation 


67.097311 

0.008413 

0.9953 

0.0125 


Regression 


Degrees 

of  Type  I  Sum 

Freedom  of  Squares  R-Square  F-Ratio  Prob  >  F 


Linear 


2  0.018495 


0.4105  130.6  0.0012 
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Quadratic 

2 

0.019983  0.4436 

141.2 

0.0011 

Crossproduct 

1 

0.006360  0.1412 

89.850 

0.0025 

Total  Regress 

5 

0.044838  0.9953 

126.7 

0.0011 

The  SAS  System 

9 

17:11  Thursday,  October 

26,  1995 

Degrees 

of  Sum  of 

Residual 

Freedom  Squares 

Mean  Square 

Total  Error 

3  0.000212 

0.000070786 

Degrees 

of 

Parameter 

Standard 

T  for  HO: 

Parameter 

Freedom 

Estimate 

Error 

Parameter^O 

Prob  >  1 T 

INTERCEPT 

1 

67.033701 

0.007201 

9308.9 

0.0000 

ALPHAL 

1 

-0.033430 

0.003348 

-9.985 

0.0021 

RATIO 

1 

-0.006629 

0.003348 

-1.980 

0.1421 

ALPHAL* ALPHAL 

1 

0.079230 

0.005135 

15.430 

0.0006 

RATIO*ALPHAL 

1 

0.039875 

0.004207 

9.479 

0.0025 

RATIO*RATIO 

1 

0.000977 

0.005135 

0.190 

0.8613 

The  SAS  System 

10 

17:11  Thursday,  October 

■  26,  1995 

Degrees 

of 

Sum  of 

Factor 

Freedom 

Squares 

Mean  Square 

F-Ratio 

Prob  >  F 

ALPHAL 

3 

0.039926 

0.013309 

188.0 

0.0007 

RATIO 

3 

0.006651 

0.002217 

31.318 

0.0092 

The  SAS  System 

11 

17:11  Thursday,  October  26,  1995 
Canonical  Analysis  of  Response  Surface 


Critical 
Factor  Value 


78 


ALPHAL 

RATIO 


0.155432 

0.220709 


Predicted  value  at  stationary  point  67.030371 


Eigenvalues 


Eigenvectors 
ALPHAL  RATIO 


0.084017  0.972366  0.233460 

-0.003810  -0.233460  0.972366 


Stationary  point  is  a  saddle  point. 

The  SAS  System  12 

17:11  Thursday,  October  26,  1995 

Estimated  Ridge  of  Minimum  Response  for  Variable  Y 


Estimated 

Standard 

Factor  Values 

Radius 

Response 

Error 

ALPHAL 

RATIO 

0 

67.040977 

0.007176 

-0.207000 

0.207000 

0.100000 

67.035841 

0.007196 

-0.110072 

0.231596 

0.200000 

67.032385 

0.007149 

-0.013603 

0.257968 

0.300000 

67.030605 

0.007027 

0.078622 

0.298760 

0.400000 

67.030087 

0.006669 

0.068286 

0.497203 

0.500000 

67.029660 

0.006317 

0.034885 

0.644598 

0.600000 

67.029172 

0.005996 

0.006015 

0.767914 

0.700000 

67.028613 

0.005734 

-0.020919 

0.881814 

0.800000 

67.027981 

0.005577 

-0.046816 

0.990799 

0.900000 

67.027274 

0.005583 

-0.072079 

1.096829 

1.000000 

67.026492 

0.005804 

-0.096925 

1.200923 

1.100000 

67.025634 

0.006271 

-0.121479 

1.303670 

1.200000 

67.024701 

0.006988 

-0.145820 

1.405439 

1.300000 

67.023692 

0.007937 

-0.170000 

1.506473 

1.400000 

67.022607 

0.009094 

-0.194056 

1.606940 

1.500000 

67.021446 

0.010434 

-0.218015 

1.706960 

The 

SAS  System 

13 

79 


17:11  Thursday,  October  26,  1995 


Estimated 

Standard 

Factor  Values 

Radius 

Response 

Error 

ALPHAL 

RATIO 

1.600000 

67.020209 

0.011936 

-0.241894 

1.806619 

1.700000 

67.018897 

0.013586 

-0.265709 

1.905986 

1.800000 

67.017508 

0.015372 

-0.289471 

2.005110 

1.900000 

67.016043 

0.017287 

-0.313188 

2.104030 

2.000000 

67.014501 

0.019323 

-0.336868 

2.202779 

2.100000 

67.012884 

0.021477 

-0.360515 

2.301381 

2.200000 

67.011191 

0.023745 

-0.384134 

2.399857 

2.300000 

67.009421 

0.026125 

-0.407729 

2.498224 

2.400000 

67.007576 

0.028614 

-0.431304 

2.596495 

2.500000 

67.005654 

0.031212 

-0.454860 

2.694683 

2.600000 

67.003656 

0.033916 

-0.478399 

2.792796 

2.700000 

67.001582 

0.036727 

-0.501924 

2.890844 

2.800000 

66.999431 

0.039643 

-0.525436 

2.988834 

2.900000 

66.997205 

0.042664 

-0.548937 

3.086771 

3.000000 

66.994902 

0.045790 

-0.572427 

3.184661 

3.100000 

66.992523 

0.049019 

-0.595908 

3.282508 

3.200000 

66.990068 

0.052351 

-0.619380 

3.380317 

3.300000 

66.987537 

0.055787 

-0.642845 

3.478091 

3.400000 

66.984929 

0.059326 

-0.666303 

3.575834 

The 

SAS  System 

14 

17 

:11  Thursday,  October  26,  1995 

Estimated 

Standard 

Factor  Values 

Radius 

Response 

Error 

ALPHAL 

RATIO 

3.500000 

66.982246 

0.062968 

-0.689754 

3.673547 

3.600000 

66.979486 

0.066712 

-0.713200 

3.771234 

3.700000 

66.976650 

0.070559 

-0.736640 

3.868896 

3.800000 

66.973738 

0.074508 

-0.760075 

3.966536 

3.900000 

66.970750 

0.078559 

-0.783505 

4.064155 

4.000000 

66.967685 

0.082713 

-0.806931 

4.161754 

4.100000 

66.964544 

0.086968 

-0.830354 

4.259336 

4.200000 

66.961327 

0.091326 

-0.853772 

4.356902 

4.300000 

66.958034 

0.095785 

-0.877188 

4.454452 

4.400000 

66.954665 

0.100346 

-0.900600 

4.551988 

4.500000 

66.951219 

0.105009 

-0.924009 

4.649510 

4.600000 

66.947698 

0.109774 

-0.947415 

4.747020 

4.700000 

4.800000 

4.900000 

5.000000 


66.944100 

66.940425 

66.936675 

66.932849 


0.114641 

0.119609 

0.124679 

0.129851 


-0.970819 

-0.994221 

-1.017620 

-1.041017 


4.844519 

4.942006 

5.039483 

5.136951 
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Appendix  D.  FORTRAN 

D.  1  Root  Mean  Squared  Error  Criterion  Program 
program  eval 

*  this  program  evaluates  the  RMS  error  criterion  from 

*  a  given  run  of  Bioplume.  It  reads  the  actual 

*  final  contaminant  concentrations  and  heads  from  the  files 

*  water991 .final  and  snapS. final  and  compares  the  values 

*  to  the  predicted  final  values  in  the  files  heads. pre 

*  and  cone. pre.  Note  that  i=column(l-20)  and  j=row  (1-30) 

■k'k'kic’kic’kic’k’k'k'k'kic'k’ic'k’k'k'k'k-k-k'k'k'k'k 

integer  i , j , junk 

real  predhead (30 , 20) ,predbenz (30, 20) , acthead (30 , 20) , 

1  aetbenz (30,20) , sumheads , sumtrans 
open (unit=l ,file='wtrfinal.bio' , status=' old' , access= 

1  'sequential' , form=' formatted' ) 
open (unit=2 , file= ' snapS .bio ' , status= ' old ' , access= 

1  'sequential' , form=' formatted' ) 
open (unit=3,file=' HEADS. BIO' ,access=' sequential' , 

1  form=' formatted' ) 

open (unit=4,file='HPLUME. BIO' ,access=' sequential' , 

1  form=' formatted' ) 
sumheads =0 
sumtrans =0 
junk=0 

************************* 

*  Fill  the  arrays  up  with  zeros 
************************* 

do  100  i=l,20 
do  200  j=l,30 

predhead (j , i) =0 .0 
predbenz ( j , i) =0 .0 
acthead (j , i) =0 . 0 
aetbenz ( j , i) =0 . 0 
200  .  continue 
100  continue 
************************* 

*  Read  the  actual  heads  and  concentrations  into  arrays 
************************* 

300  do  400  j=l,30 

read  (1,*)  (acthead( j , i) ,i=l,20) 
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read  (2,*)  (actbenz { j , i) , i=l, 20) 

400  continue 

•k'k'k'k'k'k'k'k'k’k’k'k'k'k'k'k'k'k'k'k'k'kit'k 

*  Read  the  predicted  heads  and  concentrations  into  arrays 

•k'k'k'k'k'k'k-k'k'k'k'k'k'k'kii'k'kick'k'k'k'k'k 

500  do  550  j=l,30 

read  (unit=3, fmt=*,err=600)  (predhead( j , i) , i=l, 20) 

550  continue 
goto  500 
600  continue 

read  (unit=4 ,  fint=* ,  err=675)  junk 
do  650  j=l,30 

read  (4,*)  (predbenz { j , i) , i=l, 20) 

650  continue 
goto  600 

*  Calculate  the  RMS  error  criterion  and  print  results 

675  do  700  i=l,20 
do  800  j=l,30 

sumheads=sumheads+ (acthead( j , i) -predhead( j , i) ) **2 
sumtrans=surntrans+ (actbenz (j ,i) -predbenz ( j , i) ) **2 
800  continue 

700  continue 

RMSflow= ( sumheads /600)**0.5 
RMStrans= (sumtrans/600) **0.5 

*  print  *, 'print  out  the  predbenz  array  (should  be  all  zeros) 

*  print  *, predbenz 

*  print  *, 'print  out  the  actbenz  array:' 

*  do  1000  i=l,20 

*  do  2000  j=l,30 

*  print  *, 'actbenz (row=' , j ,' ,col=' , i, ')=' ,actbenz (j , i) 

2000  continue 

1000  continue 

*  print  *, 'sumheads=' , sumheads 

*  print  *, 'sumtrans=' ,suratrans 

*  do  3000  j=l,30 

*  do  3100  i=l,20 

*  print  *,acthead( j , i) , '  vs  '  ,predhead( j , i) 

3100  continue 

3000  continue 

*  do  4000  j=l,30 
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* 

* 

4100 

4000 


do  4100  i=l,20 

print  *,actbenz(j,i) , '  vs  ' , predbenz ( j , i ) 
continue 
continue 

print  * ,  '  RMSflow= ' ,  RMSflow 
print  * , ' RMStrans= ' , RMStrans 
print  *,'done' 
end 
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